--- /dev/null
+#@ MODIF Graph Utilitai DATE 24/05/2005 AUTEUR MCOURTOI M.COURTOIS
+# -*- coding: iso-8859-1 -*-
+# CONFIGURATION MANAGEMENT OF EDF VERSION
+# ======================================================================
+# COPYRIGHT (C) 1991 - 2004 EDF R&D WWW.CODE-ASTER.ORG
+# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
+# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
+# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
+# (AT YOUR OPTION) ANY LATER VERSION.
+#
+# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
+# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
+# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
+# GENERAL PUBLIC LICENSE FOR MORE DETAILS.
+#
+# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
+# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
+# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
+# ======================================================================
+
+# RESPONSABLE MCOURTOI M.COURTOIS
+
+import sys
+import os
+import os.path
+import string
+import re
+import types
+import time
+import Numeric
+
+# try/except pour utiliser hors aster
+try:
+ import aster
+except ImportError:
+ class fake_aster:
+ def repout(self): return '/opt/aster/outils'
+ aster=fake_aster()
+
+try:
+ from Utilitai.Utmess import UTMESS
+except ImportError:
+ def UTMESS(code,sprg,texte):
+ fmt='\n <%s> <%s> %s\n\n'
+ print fmt % (code,sprg,texte)
+
+if not sys.modules.has_key('Table'):
+ try:
+ from Utilitai import Table
+ except ImportError:
+ import Table
+
+
+# ------------------------------------------------------------------------------
+class Graph:
+ """Cette classe définit l'objet Graph pour Code_Aster.
+
+ Important : Utiliser les méthodes dédiées à la manipulation des données
+ (AjoutCourbe, ...) car elles tiennent à jour les attributs "privés"
+ relatifs aux données : NbCourbe, les extrema...
+
+ Attributs :
+ - Données de chaque courbe :
+ .Valeurs : liste des valeurs de chaque courbe, pour chaque courbe :
+ (paramètres, parties réelles [, parties imaginaires])
+ .Legendes : liste des noms de chaque courbe
+ .Labels : liste des noms des colonnes de chaque courbe
+ .Styles : liste des infices de styles de ligne
+ .Couleurs : liste des indices de couleurs
+ .Marqueurs : liste des indices de symboles/marqueurs
+ .FreqMarq : liste des fréquences des marqueurs
+ .Tri : liste du tri à effectuer sur les données ('N', 'X', 'Y',
+ 'XY' ou 'YX')
+ Pour Lignes, Couleurs, Marqueurs, FreqMarq, -1 signifie valeur par défaut
+ du traceur.
+
+ - Propriétés :
+ .Titre : titre du graphique
+ .SousTitre : sous-titre (appelé commentaire dans agraf)
+ - Axes :
+ .Min_X, .Max_X, .Min_Y, .Max_Y : bornes du tracé (une méthode permet de
+ renseigner automatiquement ces valeurs avec les extréma globaux)
+ .Legende_X, .Legende_Y : légende des axes
+ .Echelle_X, .Echelle_Y : type d'échelle (LIN, LOG)
+ .Grille_X, .Grille_Y : paramètre de la grille (pas ou fréquence au choix
+ de l'utilisateur en fonction du traceur qu'il veut utiliser)
+
+ Attributs privés (modifiés uniquement par les méthodes de la classe) :
+ .NbCourbe : nombre de courbes
+ .BBXmin, BBXmax, BBYmin, BBYmax : extrema globaux (bounding box)
+ .LastTraceArgs, LastTraceFormat : données utilisées lors du dernier tracé
+ """
+# ------------------------------------------------------------------------------
+ def __init__(self):
+ """Construction + valeurs par défaut des attributs
+ """
+ self.Valeurs = []
+ self.Legendes = []
+ self.Labels = []
+ self.Styles = []
+ self.Couleurs = []
+ self.Marqueurs = []
+ self.FreqMarq = []
+ self.Tri = []
+ self.Titre = ''
+ self.SousTitre = ''
+ self.Min_X = 1.e+99
+ self.Max_X = -1.e+99
+ self.Min_Y = 1.e+99
+ self.MinP_X = 1.e+99 # minimum > 0 pour les échelles LOG
+ self.MinP_Y = 1.e+99
+ self.Max_Y = -1.e+99
+ self.Legende_X = ''
+ self.Legende_Y = ''
+ self.Echelle_X = 'LIN'
+ self.Echelle_Y = 'LIN'
+ self.Grille_X = -1
+ self.Grille_Y = -1
+ # attributs que l'utilisateur ne doit pas modifier
+ self.NbCourbe = len(self.Valeurs)
+ self.BBXmin = self.Min_X
+ self.BBXmax = self.Max_X
+ self.BBYmin = self.Min_Y
+ self.BBYmax = self.Max_Y
+ # pour conserver les paramètres du dernier tracé
+ self.LastTraceArgs = {}
+ self.LastTraceFormat = ''
+ return
+# ------------------------------------------------------------------------------
+ def SetExtrema(self,marge=0., x0=None, x1=None, y0=None, y1=None):
+ """Remplit les limites du tracé (Min/Max_X/Y) avec les valeurs de la
+ bounding box +/- avec une 'marge'*(Max-Min)/2.
+ x0,x1,y0,y1 permettent de modifier la bb.
+ """
+ if x0<>None: self.BBXmin=min([self.BBXmin, x0])
+ if x1<>None: self.BBXmax=max([self.BBXmax, x1])
+ if y0<>None: self.BBYmin=min([self.BBYmin, y0])
+ if y1<>None: self.BBYmax=max([self.BBYmax, y1])
+
+ dx=max(self.BBXmax-self.BBXmin,0.01*self.BBXmax)
+ self.Min_X = self.BBXmin - marge*dx/2.
+ self.Max_X = self.BBXmax + marge*dx/2.
+ dy=max(self.BBYmax-self.BBYmin,0.01*self.BBYmax)
+ self.Min_Y = self.BBYmin - marge*dy/2.
+ self.Max_Y = self.BBYmax + marge*dy/2.
+ return
+# ------------------------------------------------------------------------------
+ def AutoBB(self,debut=-1):
+ """Met à jour automatiquement la "bounding box"
+ (extrema toutes courbes confondues)
+ Appelé par les méthodes de manipulation des données
+ """
+ if debut == -1:
+ debut=self.NbCourbe-1
+ if debut == 0:
+ X0 = 1.e+99
+ X1 = -1.e+99
+ Y0 = 1.e+99
+ Y1 = -1.e+99
+ else:
+ X0 = self.BBXmin
+ X1 = self.BBXmax
+ Y0 = self.BBYmin
+ Y1 = self.BBYmax
+
+ for i in range(debut,self.NbCourbe):
+ X0 = min([X0,]+list(self.Valeurs[i][0]))
+ X1 = max([X1,]+list(self.Valeurs[i][0]))
+ self.MinP_X = min([self.MinP_X,]+[x for x \
+ in list(self.Valeurs[i][0]) if x>0])
+ for ny in range(1,len(self.Valeurs[i])):
+ Y0 = min([Y0,]+list(self.Valeurs[i][ny]))
+ Y1 = max([Y1,]+list(self.Valeurs[i][ny]))
+ self.MinP_Y = min([self.MinP_Y,]+[y for y \
+ in list(self.Valeurs[i][0]) if y>0])
+ self.BBXmin = X0
+ self.BBXmax = X1
+ self.BBYmin = Y0
+ self.BBYmax = Y1
+ return
+# ------------------------------------------------------------------------------
+ def AjoutCourbe(self,Val,Lab,Leg='',Sty=-1,Coul=-1,Marq=-1,FreqM=-1,Tri='N'):
+ """Ajoute une courbe dans les données
+ Val : liste de 2 listes (ou 3 si complexe) : abs, ord[, imag]
+ Leg : une chaine
+ Lab : liste de 2 chaines (ou 3 si complexe)
+ Sty : un entier
+ Coul : un entier
+ Marq : un entier
+ FreqM : un entier
+ Tri : chaine de caractères : N, X, Y, XY ou YX
+ Met à jour les attributs : NbCourbe, BBXmin/Xmax/Ymin/Ymax
+ """
+ nbc = len(Val) # nombre de colonnes : 2 ou 3
+
+ # verifications : "if not (conditions requises)"
+ if not ( 2 <= nbc <= 3 and \
+ type(Val[0]) in (types.ListType, types.TupleType) and \
+ type(Val[1]) in (types.ListType, types.TupleType) and \
+ (nbc==2 or type(Val[2]) in (types.ListType, types.TupleType)) and \
+ len(Val[0]) == len(Val[1]) and (nbc==2 or len(Val[0]) == len(Val[2])) ):
+ UTMESS('S','Graph','"Val" doit etre une liste de 2 ou 3 listes de rééls de meme longueur')
+
+ if len(Lab) <> nbc:
+ UTMESS('S','Graph','"Lab" doit etre une liste de 2 ou 3 chaines')
+
+ # ajout dans les données
+ self.Legendes.append(str(Leg))
+ self.Labels.append([str(L) for L in Lab])
+ self.Valeurs.append(Val)
+ self.Styles.append(Sty)
+ self.Couleurs.append(Coul)
+ self.Marqueurs.append(Marq)
+ self.FreqMarq.append(FreqM)
+ self.Tri.append(Tri)
+
+ self.NbCourbe = self.NbCourbe + 1
+ self.AutoBB()
+ return
+# ------------------------------------------------------------------------------
+ def Courbe(self,n):
+ """Permet de récupérer les données de la courbe d'indice n sous forme
+ d'un dictionnaire.
+ """
+ dico={
+ 'Leg' : self.Legendes[n], # légende de la courbe
+ 'LabAbs' : self.Labels[n][0], # labels des abscisses
+ 'LabOrd' : [self.Labels[n][1],], # labels des ordonnées
+ 'NbCol' : len(self.Valeurs[n]), # nombre de colonnes
+ 'NbPts' : len(self.Valeurs[n][0]), # nombre de points
+ 'Abs' : self.Valeurs[n][0], # liste des abscisses
+ 'Ord' : [self.Valeurs[n][1],], # liste des ordonnées
+ 'Sty' : self.Styles[n], # style de la ligne
+ 'Coul' : self.Couleurs[n], # couleur
+ 'Marq' : self.Marqueurs[n], # marqueur
+ 'FreqM' : self.FreqMarq[n], # fréquence du marqueur
+ 'Tri' : self.Tri[n], # ordre de tri des données
+ }
+ if(dico['NbCol'] == 3):
+ dico['LabOrd'].append(self.Labels[n][2]) # labels de la partie imaginaire
+ dico['Ord'].append(self.Valeurs[n][2]) # liste des ordonnées partie imaginaire
+ return dico
+# ------------------------------------------------------------------------------
+ def Trace(self,FICHIER=None,FORMAT=None,dform=None,**opts):
+ """Tracé du Graph selon le format spécifié.
+ FICHIER : nom du(des) fichier(s). Si None, on dirige vers stdout
+ dform : dictionnaire de formats d'impression (format des réels,
+ commentaires, saut de ligne...)
+ opts : voir TraceGraph.
+ """
+ para={
+ 'TABLEAU' : { 'mode' : 'a', 'driver' : TraceTableau, },
+ 'XMGRACE' : { 'mode' : 'a', 'driver' : TraceXmgrace, },
+ 'AGRAF' : { 'mode' : 'a', 'driver' : TraceAgraf, },
+ }
+ kargs={}
+ if self.LastTraceArgs=={}:
+ kargs['FICHIER']=FICHIER
+ kargs['dform']=dform
+ kargs['opts']=opts
+ else:
+ kargs=self.LastTraceArgs.copy()
+ if FORMAT==None:
+ FORMAT=self.LastTraceFormat
+ if FICHIER<>None:
+ kargs['FICHIER']=FICHIER
+ if dform<>None:
+ kargs['dform']=dform
+ if opts<>{}:
+ kargs['opts']=opts
+ if not FORMAT in para.keys():
+ print ' <A> <Objet Graph> Format inconnu : %s' % FORMAT
+ else:
+ kargs['fmod']=para[FORMAT]['mode']
+ self.LastTraceArgs = kargs.copy()
+ self.LastTraceFormat = FORMAT
+ # call the associated driver
+ para[FORMAT]['driver'](self,**kargs)
+# ------------------------------------------------------------------------------
+ def __repr__(self):
+ """Affichage du contenu d'un Graph"""
+ srep=''
+ for attr in ['NbCourbe','Legendes','Labels','Valeurs','Min_X','Max_X','Min_Y','Max_Y','BBXmax','BBXmin','BBYmax','BBYmin','Legende_X','Legende_Y','Echelle_X','Echelle_Y','Grille_X','Grille_Y','Tri']:
+ srep=srep + '%-10s : %s\n' % (attr,str(getattr(self,attr)))
+ return srep
+
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+class TraceGraph:
+ """
+ Cette classe définit le tracé d'un objet Graph dans un fichier.
+
+ Attributs :
+ .NomFich : liste de noms de fichier de sortie
+
+ Attributs privés (modifiés uniquement par les méthodes de la classe) :
+ .Fich : liste des objets 'fichier'
+ .Graph : objet Graph que l'on veut tracer
+ .DicForm : dictionnaire des formats de base (séparateur, format des réels...)
+
+ Les méthodes Entete, DescrCourbe, Trace (définition de l'entete, partie descriptive
+ d'une courbe, méthode de tracé/impression) sont définies dans une classe dérivée.
+ """
+# ------------------------------------------------------------------------------
+ def __init__(self,graph,FICHIER,fmod='w',dform=None,opts={}):
+ """Construction, ouverture du fichier, surcharge éventuelle du formatage
+ (dform), mode d'ouverture du fichier (fmod).
+ opts : dictionnaire dont les valeurs seront affectées comme attributs
+ de l'objet (A utiliser pour les propriétés spécifiques
+ à un format, exemple 'PILOTE' pour Xmgrace).
+ """
+ # attributs optionnels (au début pour éviter un écrasement maladroit !)
+ for k,v in opts.items():
+ setattr(self,k,v)
+
+ # Ouverture du(des) fichier(s)
+ self.NomFich=[]
+ if type(FICHIER) is types.StringType:
+ self.NomFich.append(FICHIER)
+ elif type(FICHIER) in (types.ListType, types.TupleType):
+ self.NomFich=FICHIER[:]
+ else:
+ # dans ce cas, on écrira sur stdout (augmenter le 2 éventuellement)
+ self.NomFich=[None]*2
+ self.Fich=[]
+ for ff in self.NomFich:
+ if ff<>None:
+ self.Fich.append(open(ff,fmod))
+ else:
+ self.Fich.append(sys.stdout)
+
+ # objet Graph sous-jacent
+ self.Graph=graph
+ # si Min/Max incohérents
+ if graph.Min_X > graph.Max_X or graph.Min_Y > graph.Max_Y:
+ graph.SetExtrema(marge=0.05)
+ if graph.Min_X < 0. and graph.Echelle_X=='LOG':
+ graph.Min_X=graph.MinP_X
+ if graph.Min_Y < 0. and graph.Echelle_Y=='LOG':
+ graph.Min_Y=graph.MinP_Y
+
+ # formats de base (identiques à ceux du module Table)
+ self.DicForm={
+ 'csep' : ' ', # séparateur
+ 'ccom' : '#', # commentaire
+ 'cdeb' : '', # début de ligne
+ 'cfin' : '\n', # fin de ligne
+ 'formK' : '%-12s', # chaines
+ 'formR' : '%12.5E', # réels
+ 'formI' : '%12d' # entiers
+ }
+ if dform<>None and type(dform)==types.DictType:
+ self.DicForm.update(dform)
+
+ # let's go
+ self.Trace()
+ return
+# ------------------------------------------------------------------------------
+ def __del__(self):
+ """Fermeture du(des) fichier(s) à la destruction"""
+ if hasattr(self,'Fich'):
+ self._FermFich()
+# ------------------------------------------------------------------------------
+ def _FermFich(self):
+ """Fermeture du(des) fichier(s)"""
+ for fp in self.Fich:
+ if fp<>sys.stdout:
+ fp.close()
+# ------------------------------------------------------------------------------
+ def _OuvrFich(self):
+ """Les fichiers sont ouverts par le constructeur. S'ils ont été fermés,
+ par un appel au Tracé, _OuvrFich ouvre de nouveau les fichiers dans le
+ meme mode"""
+ n=len(self.NomFich)
+ for i in range(n):
+ if self.Fich[i].closed:
+ self.Fich[i]=open(self.NomFich[i],self.Fich[i].mode)
+
+# ------------------------------------------------------------------------------
+ def Entete(self):
+ """Retourne l'entete"""
+ raise StandardError, "Cette méthode doit etre définie par la classe fille."
+# ------------------------------------------------------------------------------
+ def DescrCourbe(self,**args):
+ """Retourne la chaine de caractères décrivant les paramètres de la courbe.
+ """
+ raise StandardError, "Cette méthode doit etre définie par la classe fille."
+# ------------------------------------------------------------------------------
+ def Trace(self):
+ """Méthode pour 'tracer' l'objet Graph dans un fichier.
+ Met en page l'entete, la description des courbes et les valeurs selon
+ le format et ferme le fichier.
+ """
+ raise StandardError, "Cette méthode doit etre définie par la classe fille."
+
+
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+class TraceTableau(TraceGraph):
+ """
+ Impression d'un objet Graph sous forme d'un tableau de colonnes,
+ on suppose que les courbes partagent la meme liste d'abscisse à 'EPSILON'
+ près, sinon on alarme.
+ """
+ EPSILON=1.e-4
+# ------------------------------------------------------------------------------
+ def Trace(self):
+ """Méthode pour 'tracer' l'objet Graph dans un fichier.
+ Met en page l'entete, la description des courbes et les valeurs selon
+ le format et ferme le fichier.
+ L'ouverture et la fermeture du fichier sont gérées par l'objet Table.
+ """
+ g=self.Graph
+ msg=[]
+ if g.NbCourbe > 0:
+ # validité des données (abscisses identiques)
+ t0=Numeric.array(g.Courbe(0)['Abs'])
+ max0=max(abs(t0))
+ for i in range(1,g.NbCourbe):
+ if g.Courbe(i)['NbPts']<>g.Courbe(0)['NbPts']:
+ msg.append(" <A> <TraceTableau> La courbe %d n'a pas le meme " \
+ "nombre de points que la 1ère." % i)
+ else:
+ ti=Numeric.array(g.Courbe(i)['Abs'])
+ if max(abs((ti-t0).flat)) > self.EPSILON*max0:
+ msg.append(" <A> <TraceTableau> Courbe %d : écart entre les "\
+ "abscisses supérieur à %9.2E" % (i+1,self.EPSILON))
+ msg.append(" Utilisez IMPR_FONCTION pour interpoler " \
+ "les valeurs sur la première liste d'abscisses.")
+ # objet Table
+ Tab=Table.Table()
+ # titre / sous-titre
+ tit=[]
+ tit.append(self.DicForm['ccom']+' '+g.Titre)
+ tit.append(self.DicForm['ccom']+' '+g.SousTitre)
+ # legendes
+ for i in range(g.NbCourbe):
+ tit.append(self.DicForm['ccom']+' Courbe '+str(i)+' '+g.Legendes[i])
+ Tab.titr=self.DicForm['cfin'].join(tit)
+ # noms des paramètres/colonnes
+ Tab.para.append(g.Labels[0][0])
+ for i in range(g.NbCourbe):
+ for lab in g.Labels[i][1:]:
+ Tab.para.append(lab)
+ # types
+ Tab.type=['R']*len(Tab.para)
+ # lignes de la Table
+ dC0=g.Courbe(0)
+ for j in range(dC0['NbPts']):
+ row={}
+ row[dC0['LabAbs']]=dC0['Abs'][j]
+ for i in range(g.NbCourbe):
+ dCi=g.Courbe(i)
+ for k in range(dCi['NbCol']-1):
+ try:
+ row[dCi['LabOrd'][k]]=dCi['Ord'][k][j]
+ except IndexError:
+ row[dCi['LabOrd'][k]]=None
+ Tab.append(row)
+ Tab.Impr(FICHIER=self.NomFich[0], FORMAT='TABLEAU')
+ # erreurs ?
+ if msg:
+ print '\n'.join(msg)
+ return
+
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+class TraceXmgrace(TraceGraph):
+ """
+ Impression d'un objet Graph au format XMGRACE.
+ Attribut supplémentaire : .PILOTE
+ """
+ PILOTE=''
+# ------------------------------------------------------------------------------
+ def Entete(self):
+ """Retourne l'entete du fichier .agr correspondant à la mise en forme
+ """
+ dic_ech={ 'LIN' : 'Normal', 'LOG' : 'Logarithmic' }
+ g=self.Graph
+ entete=[]
+ entete.append("""
+# Grace project file
+#
+@version 50100
+@page size 842, 595
+@page scroll 5%
+@page inout 5%
+@link page off
+@map font 0 to "Times-Roman", "Times-Roman"
+@map font 1 to "Times-Italic", "Times-Italic"
+@map font 2 to "Times-Bold", "Times-Bold"
+@map font 3 to "Times-BoldItalic", "Times-BoldItalic"
+@map font 4 to "Helvetica", "Helvetica"
+@map font 5 to "Helvetica-Oblique", "Helvetica-Oblique"
+@map font 6 to "Helvetica-Bold", "Helvetica-Bold"
+@map font 7 to "Helvetica-BoldOblique", "Helvetica-BoldOblique"
+@map font 8 to "Courier", "Courier"
+@map font 9 to "Courier-Oblique", "Courier-Oblique"
+@map font 10 to "Courier-Bold", "Courier-Bold"
+@map font 11 to "Courier-BoldOblique", "Courier-BoldOblique"
+@map font 12 to "Symbol", "Symbol"
+@map font 13 to "ZapfDingbats", "ZapfDingbats"
+@map color 0 to (255, 255, 255), "white"
+@map color 1 to (0, 0, 0), "black"
+@map color 2 to (255, 0, 0), "red"
+@map color 3 to (0, 255, 0), "green"
+@map color 4 to (0, 0, 255), "blue"
+@map color 5 to (255, 255, 0), "yellow"
+@map color 6 to (188, 143, 143), "brown"
+@map color 7 to (220, 220, 220), "grey"
+@map color 8 to (148, 0, 211), "violet"
+@map color 9 to (0, 255, 255), "cyan"
+@map color 10 to (255, 0, 255), "magenta"
+@map color 11 to (255, 165, 0), "orange"
+@map color 12 to (114, 33, 188), "indigo"
+@map color 13 to (103, 7, 72), "maroon"
+@map color 14 to (64, 224, 208), "turquoise"
+@map color 15 to (0, 139, 0), "green4"
+@reference date 0
+@date wrap off
+@date wrap year 1950
+@timestamp off
+@default linewidth 1.0
+@default linestyle 1
+@default color 1
+@default pattern 1
+@default font 0
+@default char size 1.000000
+@default symbol size 1.000000
+@default sformat "%.8g"
+@background color 0
+@page background fill on
+@r0 off
+@link r0 to g0
+@r0 type above
+@r0 linestyle 1
+@r0 linewidth 1.0
+@r0 color 1
+@r0 line 0, 0, 0, 0
+@r1 off
+@link r1 to g0
+@r1 type above
+@r1 linestyle 1
+@r1 linewidth 1.0
+@r1 color 1
+@r1 line 0, 0, 0, 0
+@r2 off
+@link r2 to g0
+@r2 type above
+@r2 linestyle 1
+@r2 linewidth 1.0
+@r2 color 1
+@r2 line 0, 0, 0, 0
+@r3 off
+@link r3 to g0
+@r3 type above
+@r3 linestyle 1
+@r3 linewidth 1.0
+@r3 color 1
+@r3 line 0, 0, 0, 0
+@r4 off
+@link r4 to g0
+@r4 type above
+@r4 linestyle 1
+@r4 linewidth 1.0
+@r4 color 1
+@r4 line 0, 0, 0, 0
+@g0 on
+@g0 hidden false
+@g0 type XY
+@g0 stacked false
+@g0 bar hgap 0.000000
+@with g0
+@ stack world 0, 0, 0, 0
+@ znorm 1
+@ view xmin 0.150000
+@ view xmax 1.150000
+@ view ymin 0.150000
+@ view ymax 0.850000
+@ title font 0
+@ title size 1.500000
+@ title color 1
+@ subtitle font 0
+@ subtitle size 1.000000
+@ subtitle color 1
+@ xaxes invert off
+@ yaxes invert off
+@ xaxis on
+@ xaxis type zero false
+@ xaxis offset 0.000000 , 0.000000
+@ xaxis bar on
+@ xaxis bar color 1
+@ xaxis bar linestyle 1
+@ xaxis bar linewidth 1.0
+@ xaxis label layout para
+@ xaxis label place auto
+@ xaxis label char size 1.000000
+@ xaxis label font 0
+@ xaxis label color 1
+@ xaxis label place normal
+@ xaxis tick on
+@ xaxis tick minor ticks 1
+@ xaxis tick default 6
+@ xaxis tick place rounded true
+@ xaxis tick in
+@ xaxis tick major size 1.000000
+@ xaxis tick major color 1
+@ xaxis tick major linewidth 1.0
+@ xaxis tick major linestyle 2
+@ xaxis tick major grid on
+@ xaxis tick minor color 1
+@ xaxis tick minor linewidth 1.0
+@ xaxis tick minor linestyle 2
+@ xaxis tick minor grid off
+@ xaxis tick minor size 0.500000
+@ xaxis ticklabel on
+@ xaxis ticklabel format general
+@ xaxis ticklabel prec 5
+@ xaxis ticklabel angle 0
+@ xaxis ticklabel skip 0
+@ xaxis ticklabel stagger 0
+@ xaxis ticklabel place normal
+@ xaxis ticklabel offset auto
+@ xaxis ticklabel offset 0.000000 , 0.010000
+@ xaxis ticklabel start type auto
+@ xaxis ticklabel start 0.000000
+@ xaxis ticklabel stop type auto
+@ xaxis ticklabel stop 0.000000
+@ xaxis ticklabel char size 0.800000
+@ xaxis ticklabel font 0
+@ xaxis ticklabel color 1
+@ xaxis ticklabel formula ""
+@ xaxis ticklabel append ""
+@ xaxis ticklabel prepend ""
+@ xaxis tick place both
+@ xaxis tick spec type none
+@ yaxis on
+@ yaxis type zero false
+@ yaxis offset 0.000000 , 0.000000
+@ yaxis bar on
+@ yaxis bar color 1
+@ yaxis bar linestyle 1
+@ yaxis bar linewidth 1.0
+@ yaxis label layout para
+@ yaxis label place auto
+@ yaxis label char size 1.000000
+@ yaxis label font 0
+@ yaxis label color 1
+@ yaxis label place normal
+@ yaxis tick on
+@ yaxis tick minor ticks 1
+@ yaxis tick default 6
+@ yaxis tick place rounded true
+@ yaxis tick in
+@ yaxis tick major size 1.000000
+@ yaxis tick major color 1
+@ yaxis tick major linewidth 1.0
+@ yaxis tick major linestyle 2
+@ yaxis tick major grid on
+@ yaxis tick minor color 1
+@ yaxis tick minor linewidth 1.0
+@ yaxis tick minor linestyle 1
+@ yaxis tick minor grid off
+@ yaxis tick minor size 0.500000
+@ yaxis ticklabel on
+@ yaxis ticklabel format general
+@ yaxis ticklabel prec 5
+@ yaxis ticklabel angle 0
+@ yaxis ticklabel skip 0
+@ yaxis ticklabel stagger 0
+@ yaxis ticklabel place normal
+@ yaxis ticklabel offset auto
+@ yaxis ticklabel offset 0.000000 , 0.010000
+@ yaxis ticklabel start type auto
+@ yaxis ticklabel start 0.000000
+@ yaxis ticklabel stop type auto
+@ yaxis ticklabel stop 0.000000
+@ yaxis ticklabel char size 0.800000
+@ yaxis ticklabel font 0
+@ yaxis ticklabel color 1
+@ yaxis ticklabel formula ""
+@ yaxis ticklabel append ""
+@ yaxis ticklabel prepend ""
+@ yaxis tick place both
+@ yaxis tick spec type none
+@ altxaxis off
+@ altyaxis off
+@ legend on
+@ legend loctype view
+@ legend 0.85, 0.8
+@ legend box color 1
+@ legend box pattern 1
+@ legend box linewidth 1.0
+@ legend box linestyle 1
+@ legend box fill color 0
+@ legend box fill pattern 1
+@ legend font 0
+@ legend char size 0.750000
+@ legend color 1
+@ legend length 4
+@ legend vgap 1
+@ legend hgap 1
+@ legend invert false
+@ frame type 0
+@ frame linestyle 1
+@ frame linewidth 1.0
+@ frame color 1
+@ frame pattern 1
+@ frame background color 0
+@ frame background pattern 0
+""")
+ entete.append('@ title "'+g.Titre+'"')
+ entete.append('@ subtitle "'+g.SousTitre+'"')
+ entete.append('@ xaxis label "'+g.Legende_X+'"')
+ entete.append('@ yaxis label "'+g.Legende_Y+'"')
+ entete.append('@ xaxes scale '+dic_ech[g.Echelle_X])
+ entete.append('@ yaxes scale '+dic_ech[g.Echelle_Y])
+ entete.append('@ xaxis tick major '+str(g.Grille_X))
+ entete.append('@ yaxis tick major '+str(g.Grille_Y))
+ entete.append('@ world xmin '+str(g.Min_X))
+ entete.append('@ world xmax '+str(g.Max_X))
+ entete.append('@ world ymin '+str(g.Min_Y))
+ entete.append('@ world ymax '+str(g.Max_Y))
+ return entete
+# ------------------------------------------------------------------------------
+ def DescrCourbe(self,**args):
+ """Retourne la chaine de caractères décrivant les paramètres de la courbe.
+ """
+ # valeurs par défaut
+ sty = str(ValCycl(args['Sty'],0,8,1))
+ color = str(ValCycl(args['Coul'],1,15,args['NumSet']+1))
+ symbol= str(ValCycl(args['Marq'],0,10,args['NumSet']))
+ freqm = str(ValCycl(args['FreqM'],0,-1,0))
+
+ sn=str(args['NumSet'])
+ descr=[]
+ descr.append(string.replace("""
+@ s0 hidden false
+@ s0 type xy
+@ s0 symbol size 1.000000
+@ s0 symbol pattern 1
+@ s0 symbol linestyle 1
+@ s0 symbol fill pattern 0
+@ s0 symbol linewidth 1.0
+@ s0 symbol char 65
+@ s0 symbol char font 0
+@ s0 line type 1
+@ s0 line linewidth 1.0
+@ s0 line pattern 1
+@ s0 baseline type 0
+@ s0 baseline off
+@ s0 dropline off
+@ s0 fill type 0
+@ s0 fill rule 0
+@ s0 fill pattern 1
+@ s0 avalue off
+@ s0 avalue type 2
+@ s0 avalue char size 1.000000
+@ s0 avalue font 0
+@ s0 avalue rot 0
+@ s0 avalue format general
+@ s0 avalue prec 3
+@ s0 avalue prepend ""
+@ s0 avalue append ""
+@ s0 avalue offset 0.000000 , 0.000000
+@ s0 errorbar on
+@ s0 errorbar place both
+@ s0 errorbar pattern 1
+@ s0 errorbar size 1.000000
+@ s0 errorbar linewidth 1.0
+@ s0 errorbar linestyle 1
+@ s0 errorbar riser linewidth 1.0
+@ s0 errorbar riser linestyle 1
+@ s0 errorbar riser clip off
+@ s0 errorbar riser clip length 0.100000
+
+@ s0 comment ""
+""",' s0 ',' s'+sn+' '))
+ descr.append('@ s'+sn+' symbol '+symbol)
+ descr.append('@ s'+sn+' symbol color '+color)
+ descr.append('@ s'+sn+' symbol skip '+freqm)
+ descr.append('@ s'+sn+' symbol fill color '+color)
+ descr.append('@ s'+sn+' line linestyle '+sty)
+ descr.append('@ s'+sn+' line color '+color)
+ descr.append('@ s'+sn+' fill color '+color)
+ descr.append('@ s'+sn+' avalue color '+color)
+ descr.append('@ s'+sn+' errorbar color '+color)
+ descr.append('@ s'+sn+' legend "'+args['Leg']+'"')
+ return descr
+# ------------------------------------------------------------------------------
+ def Trace(self):
+ """Méthode pour 'tracer' l'objet Graph dans un fichier.
+ Met en page l'entete, la description des courbes et les valeurs selon
+ le format et ferme le fichier.
+ """
+ g=self.Graph
+ if self.PILOTE=='INTERACTIF':
+ self.NomFich[0]='Trace_'+time.strftime('%y%m%d%H%M%S',time.localtime())+'.dat'
+ self.Fich[0]=open(self.NomFich[0],'w')
+ # initialise le graph
+ self._FermFich()
+ nbsets, x0, x1, y0, y1 = IniGrace(self.NomFich[0])
+ NumSetIni = nbsets+1
+ g.SetExtrema(0.05, x0, x1, y0, y1)
+ # si Min/Max incohérents
+ if g.Min_X < 0. and g.Echelle_X=='LOG':
+ g.Min_X=g.MinP_X
+ if g.Min_Y < 0. and g.Echelle_Y=='LOG':
+ g.Min_Y=g.MinP_Y
+
+ self._OuvrFich()
+ fich=self.Fich[0]
+ if g.NbCourbe < 1:
+ self._FermFich()
+ return
+ # cohérence des valeurs par défaut
+ if g.Grille_X<0 or g.Grille_Y<0:
+ deltaX=g.Max_X-g.Min_X
+ deltaY=g.Max_Y-g.Min_Y
+ g.Grille_X=deltaX/5.
+ g.Grille_Y=deltaY/5.
+ if deltaX>4:
+ g.Grille_X=int(round(g.Grille_X))
+ if deltaY>4:
+ g.Grille_Y=int(round(g.Grille_Y))
+ # entete
+ fich.write('\n'.join(self.Entete()))
+ fich.write('\n')
+ # valeurs
+ it=-1
+ for i in range(g.NbCourbe):
+ dCi=g.Courbe(i)
+ for k in range(dCi['NbCol']-1):
+ it=it+1
+ dCi['NumSet'] = NumSetIni + it
+ fich.write('\n'.join(self.DescrCourbe(**dCi)))
+ fich.write('\n')
+ # partie données (.dat)
+ lig=[]
+ it=-1
+ for i in range(g.NbCourbe):
+ dCi=g.Courbe(i)
+ for k in range(dCi['NbCol']-1):
+ it=it+1
+ lig.append('@target g0.s%d' % (NumSetIni + it))
+ lig.append('@type xy')
+ listX, listY = Tri(g.Tri, lx=dCi['Abs'], ly=dCi['Ord'][k])
+ for j in range(dCi['NbPts']):
+ svX=self.DicForm['formR'] % listX[j]
+ svY=self.DicForm['formR'] % listY[j]
+ lig.append(self.DicForm['formR'] % listX[j] + \
+ ' ' + self.DicForm['formR'] % listY[j])
+ lig.append('&')
+ fich.write('\n'.join(lig))
+ fich.write('\n')
+ self._FermFich()
+
+ # Production du fichier postscript, jpeg ou lancement interactif
+ pilo=self.PILOTE
+ if self.PILOTE<>'':
+ xmgr=os.path.join(aster.repout(),'xmgrace')
+ nfhard=self.NomFich[0]+'.hardcopy'
+ # nom exact du pilote
+ if pilo=='POSTSCRIPT':
+ pilo='PostScript'
+ elif pilo=='INTERACTIF':
+ pilo='X11'
+ # ligne de commande
+ if pilo=='X11':
+ lcmde=xmgr+' '+self.NomFich[0]
+ if not os.environ.has_key('DISPLAY') or os.environ['DISPLAY']=='':
+ os.environ['DISPLAY']=':0.0'
+ UTMESS('A','TraceXmgrace','Variable DISPLAY non définie')
+ UTMESS('I','TraceXmgrace','on fixe le DISPLAY à %s' % os.environ['DISPLAY'])
+ else:
+ if os.path.exists(os.path.join(aster.repout(),'gracebat')):
+ xmgr=os.path.join(aster.repout(),'gracebat')
+ lcmde=xmgr+' -hdevice '+pilo+' -hardcopy -printfile '+nfhard+' '+self.NomFich[0]
+ # appel xmgrace
+ UTMESS('I','TraceXmgrace','Lancement de : '+lcmde)
+ if not os.path.exists(xmgr):
+ UTMESS('S','TraceXmgrace','Fichier inexistant : '+xmgr)
+ iret=os.system(lcmde)
+ if iret==0 or os.path.exists(nfhard):
+ if pilo not in ['','X11']:
+ os.remove(self.NomFich[0]) # necessaire sous windows
+ os.rename(nfhard,self.NomFich[0])
+ else:
+ UTMESS('A','TraceXmgrace',"Erreur lors de l'utilisation du filtre "+pilo+"\nLe fichier retourné est le fichier '.agr'")
+ # menage
+ if self.PILOTE=='INTERACTIF':
+ os.remove(self.NomFich[0])
+ return
+
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+class TraceAgraf(TraceGraph):
+ """
+ Impression d'un objet Graph au format AGRAF.
+ """
+# ------------------------------------------------------------------------------
+ def Entete(self):
+ """Retourne l'entete des directives Agraf"""
+ dic_ech={ 'LIN' : '0', 'LOG' : '1' }
+ g=self.Graph
+ entete=[]
+ entete.append("""
+ASPECT_GRAPHIQUE:
+ En-tete :Departement Analyses Mecaniques et Acoustique
+ Aspect :0
+ Nombre de vues :1
+ Cesure commentaire :40
+ MinMax :0
+ Fonte Titre :%helvetica-14
+ Fonte Axes :%courier-12
+ Fonte Autre :%times-12
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 0 0 0
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 65535 0 0
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 11822 35723 22359
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 0 0 65535
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 65535 0 65535
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 0 65535 65535
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 0 65535 0
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 41120 21074 11565
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 65535 42405 0
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 41120 8224 61680
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 65535 65535 0
+
+ DEFAUT_COURBE:
+ Couleur (rvb) : 53970 46260 35980
+
+GRAPHIQUE:
+""")
+ if g.Titre=='':
+ g.Titre='GRAPHIQUE CODE_ASTER'
+ entete.append('Titre :'+g.Titre+'\n')
+ if g.SousTitre<>'':
+ entete.append('Commentaire :'+g.SousTitre+'\n')
+ entete.append('Frequence Grille X :'+str(int(g.Grille_X))+'\n')
+ entete.append('Frequence Grille Y :'+str(int(g.Grille_Y))+'\n')
+ entete.append('Echelle X :'+dic_ech[g.Echelle_X]+'\n')
+ entete.append('Echelle Y :'+dic_ech[g.Echelle_Y]+'\n')
+ if g.Legende_X<>'':
+ entete.append('Legende X :'+g.Legende_X+'\n')
+ if g.Legende_Y<>'':
+ entete.append('Legende Y :'+g.Legende_Y+'\n')
+ entete.append('Min X : '+str(g.Min_X)+'\n')
+ entete.append('Max X : '+str(g.Max_X)+'\n')
+ entete.append('Min Y : '+str(g.Min_Y)+'\n')
+ entete.append('Max Y : '+str(g.Max_Y)+'\n')
+
+ return entete
+# ------------------------------------------------------------------------------
+ def DescrCourbe(self,**args):
+ """Retourne la chaine de caractères décrivant les paramètres de la courbe.
+ """
+ # valeurs par défaut
+ sty = str(ValCycl(args['Sty'],0,2,0))
+ color = str(ValCycl(args['Coul'],0,12,args['NumSet']))
+ symbol= str(ValCycl(args['Marq'],0,12,args['NumSet']))
+ freqm = str(ValCycl(args['FreqM'],0,-1,0))
+
+ descr=[]
+ descr.append(' COURBE:\n')
+ descr.append(' Trait :'+sty+'\n')
+ descr.append(' Couleur :'+color+'\n')
+ descr.append(' Marqueur :'+symbol+'\n')
+ descr.append(' Frequence Marqueur :'+freqm+'\n')
+ if args['Leg']<>'':
+ descr.append(' Legende :'+args['Leg']+'\n')
+ descr.append(' Tri :'+args['Tri']+'\n')
+ descr.append(' Abscisses : [ '+str(args['Bloc'])+', '+str(args['ColX'])+']\n')
+ descr.append(' Ordonnees : [ '+str(args['Bloc'])+', '+str(args['ColY'])+']\n')
+ return descr
+# ------------------------------------------------------------------------------
+ def Trace(self):
+ """Méthode pour 'tracer' l'objet Graph dans un fichier.
+ Met en page l'entete, la description des courbes et les valeurs selon
+ le format et ferme le fichier.
+ """
+ self._OuvrFich()
+ fdogr=self.Fich[0]
+ fdigr=self.Fich[1]
+ g=self.Graph
+ if g.NbCourbe > 0:
+ # cohérence des valeurs par défaut
+ if g.Grille_X<0 or g.Grille_Y<0:
+ g.Grille_X=0
+ g.Grille_Y=0
+ # entete
+ for lig in self.Entete():
+ fdigr.write(lig)
+ # valeurs
+ for i in range(g.NbCourbe):
+ dCi=g.Courbe(i)
+ dCi['NumSet']=i
+ # partie directives (.digr)
+ for k in range(dCi['NbCol']-1):
+ dCi['Bloc']=i+1
+ dCi['ColX']=1
+ dCi['ColY']=k+2
+ for lig in self.DescrCourbe(**dCi):
+ fdigr.write(lig)
+ # partie données (.dogr)
+ if dCi['Leg']<>'':
+ leg=dCi['Leg']
+ else:
+ leg='COURBE_'+str(i)
+ fdogr.write('#NOM DE LA FONCTION: '+leg+'\n')
+ for j in range(dCi['NbPts']):
+ for k in range(dCi['NbCol']):
+ sv=self.DicForm['formR'] % g.Valeurs[i][k][j]
+ fdogr.write(' '+sv)
+ fdogr.write('\n')
+ fdogr.write('\n')
+ self._FermFich()
+
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+def ValCycl(val,vmin,vmax,vdef):
+ """
+ Retourne une valeur entre vmin et vmax (bornes incluses) :
+ - si val<vmin, on utilise val=vdef,
+ - si val>vmax, on cycle tel que val=vmax+1 retourne vmin, etc.
+ - si vmax<vmin, il n'y a pas de max
+ """
+ v = val
+ if v < vmin:
+ v = vdef
+ if vmax < vmin:
+ return v
+ else:
+ return (((v-vmin) % (vmax+1-vmin))+vmin)
+
+# ------------------------------------------------------------------------------
+def Tri(tri, lx, ly):
+ """Retourne les listes triées selon la valeur de tri ('X','Y','XY','YX').
+ """
+ dNumCol={ 'X' : 0, 'Y' : 1 }
+ tab=Numeric.array((lx,ly),Numeric.Float64)
+ tab=Numeric.transpose(tab)
+ li=range(len(tri))
+ li.reverse()
+ for i in li:
+ if tri[-i] in dNumCol.keys():
+ icol=dNumCol[tri[-i]]
+ tab = Numeric.take(tab, Numeric.argsort(tab[:,icol]))
+ return [ tab[:,0].tolist(), tab[:,1].tolist() ]
+
+# ------------------------------------------------------------------------------
+def AjoutParaCourbe(dCourbe, args):
+ """Ajoute les arguments fournis par l'utilisateur (args) dans le dictionnaire
+ décrivant la courbe (dCourbe).
+ """
+ # correspondance : mot-clé Aster / clé du dico de l'objet Graph
+ keys={
+ 'LEGENDE' : 'Leg',
+ 'STYLE' : 'Sty',
+ 'COULEUR' : 'Coul',
+ 'MARQUEUR' : 'Marq',
+ 'FREQ_MARQUEUR' : 'FreqM',
+ 'TRI' : 'Tri',
+ }
+ for mc, key in keys.items():
+ if args.has_key(mc):
+ dCourbe[key]=args[mc]
+
+# ------------------------------------------------------------------------------
+def IniGrace(fich):
+ """Retourne le numéro de la dernière courbe d'un fichier xmgrace (sinon 0).
+ """
+ ns=0
+ x0=None
+ x1=None
+ y0=None
+ y1=None
+ if os.path.exists(fich) and os.stat(fich).st_size<>0:
+ os.rename(fich, fich+'.prev')
+ fpre=open(fich+'.prev', 'r')
+ fnew=open(fich, 'w')
+ for line in fpre:
+ ikeep=True
+ mat=re.search('@target g[0-9]+\.s([0-9]+)', line)
+ if mat<>None and int(mat.group(1))>ns:
+ ns=int(mat.group(1))
+ mat=re.search('@[ ]+world[ ]+xmin[ ]+([\-\+\.0-9eEdD]+)', line)
+ if mat<>None:
+ try:
+ x0=float(mat.group(1))
+ ikeep=False
+ except ValueError:
+ pass
+ mat=re.search('@[ ]+world[ ]+xmax[ ]+([\-\+\.0-9eEdD]+)', line)
+ if mat<>None:
+ try:
+ x1=float(mat.group(1))
+ ikeep=False
+ except ValueError:
+ pass
+ mat=re.search('@[ ]+world[ ]+ymin[ ]+([\-\+\.0-9eEdD]+)', line)
+ if mat<>None:
+ try:
+ y0=float(mat.group(1))
+ ikeep=False
+ except ValueError:
+ pass
+ mat=re.search('@[ ]+world[ ]+ymax[ ]+([\-\+\.0-9eEdD]+)', line)
+ if mat<>None:
+ try:
+ y1=float(mat.group(1))
+ ikeep=False
+ except ValueError:
+ pass
+ if ikeep:
+ fnew.write(line)
+ fpre.close()
+ fnew.close()
+ print """
+ <I> Informations sur le fichier '%s' :
+ Nombre de courbes : %3d
+ Bornes des abscisses : [ %13.6G , %13.6G ]
+ Bornes des ordonnées : [ %13.6G , %13.6G ]
+""" % (fich, ns, x0, x1, y0, y1)
+ return ns, x0, x1, y0, y1
--- /dev/null
+#@ MODIF Table Utilitai DATE 17/05/2005 AUTEUR DURAND C.DURAND
+# -*- coding: iso-8859-1 -*-
+# CONFIGURATION MANAGEMENT OF EDF VERSION
+# ======================================================================
+# COPYRIGHT (C) 1991 - 2004 EDF R&D WWW.CODE-ASTER.ORG
+# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
+# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
+# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
+# (AT YOUR OPTION) ANY LATER VERSION.
+#
+# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
+# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
+# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
+# GENERAL PUBLIC LICENSE FOR MORE DETAILS.
+#
+# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
+# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
+# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
+# ======================================================================
+
+# RESPONSABLE MCOURTOI M.COURTOIS
+
+import sys
+import string
+import re
+
+from types import *
+EnumTypes=(ListType, TupleType)
+NumberTypes=(IntType, LongType, FloatType, ComplexType)
+
+# try/except pour utiliser hors aster
+try:
+ from Utilitai.Utmess import UTMESS
+except ImportError:
+ def UTMESS(code,sprg,texte):
+ fmt='\n <%s> <%s> %s\n\n'
+ print fmt % (code,sprg,texte)
+
+if not sys.modules.has_key('Graph'):
+ try:
+ from Utilitai import Graph
+ except ImportError:
+ import Graph
+
+# formats de base (identiques à ceux du module Graph)
+DicForm={
+ 'csep' : ' ', # séparateur
+ 'ccom' : '#', # commentaire
+ 'cdeb' : '', # début de ligne
+ 'cfin' : '\n', # fin de ligne
+ 'formK' : '%-8s', # chaines
+ 'formR' : '%12.5E', # réels
+ 'formI' : '%8d' # entiers
+}
+
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+class TableBase(object):
+ """Classe pour partager les méthodes d'impression entre Table et Colonne
+ (c'est surtout utile pour vérifier que l'extraction et les filtres sur les
+ colonnes sont corrects).
+ """
+ def __repr__(self):
+ return self.ReprTable()
+ def Croise(self,**kargs):
+ raise StandardError, 'Must be defined in a derived class'
+
+# ------------------------------------------------------------------------------
+ def Impr(self,FICHIER=None,FORMAT='TABLEAU',dform=None,**opts):
+ """Impresssion de la Table selon le format spécifié.
+ FICHIER : nom du(des) fichier(s). Si None, on dirige vers stdout
+ dform : dictionnaire de formats d'impression (format des réels,
+ commentaires, saut de ligne...)
+ opts : selon FORMAT.
+ """
+ para={
+ 'TABLEAU' : { 'mode' : 'a', 'driver' : self.ImprTableau, },
+ 'ASTER' : { 'mode' : 'a', 'driver' : self.ImprTableau, },
+ 'XMGRACE' : { 'mode' : 'a', 'driver' : self.ImprGraph, },
+ 'AGRAF' : { 'mode' : 'a', 'driver' : self.ImprTableau, },
+ 'TABLEAU_CROISE' : { 'mode' : 'a', 'driver' : self.ImprTabCroise, },
+ }
+ kargs={
+ 'FICHIER' : FICHIER,
+ 'FORMAT' : FORMAT,
+ 'dform' : DicForm.copy(),
+ 'mode' : para[FORMAT]['mode'],
+ }
+ if dform<>None and type(dform)==DictType:
+ kargs['dform'].update(dform)
+ # ajout des options
+ kargs.update(opts)
+
+ if not kargs.get('PAGINATION'):
+ # call the associated driver
+ para[FORMAT]['driver'](**kargs)
+
+ else:
+ if not type(kargs['PAGINATION']) in EnumTypes:
+ ppag=[kargs['PAGINATION'],]
+ else:
+ ppag=list(kargs['PAGINATION'])
+ del kargs['PAGINATION']
+ npag=len(ppag)
+ # paramètres hors ceux de la pagination
+ lkeep=[p for p in self.para if ppag.count(p)==0]
+ # création des listes des valeurs distinctes
+ lvd=[]
+ for p in ppag:
+ lvp=getattr(self,p).values()
+ lvn=[]
+ for it in lvp:
+ if it<>None and lvn.count(it)==0:
+ lvn.append(it)
+ lvn.sort()
+ lvd.append(lvn)
+ # création des n-uplets
+ s = '[['+','.join(['x'+str(i) for i in range(npag)])+'] '
+ s+= ' '.join(['for x'+str(i)+' in lvd['+str(i)+']' for i in range(npag)])+']'
+ try:
+ lnup=eval(s)
+ except SyntaxError, s:
+ UTMESS('F','Table','Erreur lors de la construction des n-uplets')
+ # pour chaque n-uplet, on imprime la sous-table
+ for nup in lnup:
+ tab=self
+ for i in range(npag):
+ tab = tab & (getattr(tab,ppag[i]) == nup[i])
+ sl=''
+ if tab.titr: sl='\n'
+ tab.titr += sl+ppag[i]+': '+str(nup[i])
+ tab[lkeep].Impr(**kargs)
+
+# ------------------------------------------------------------------------------
+ def ImprTableau(self,**kargs):
+ """Impression au format TABLEAU ou ASTER
+ """
+ # fichier ou stdout
+ if kargs.get('FICHIER')<>None:
+ f=open(kargs['FICHIER'],kargs['mode'])
+ else:
+ f=sys.stdout
+ # ecriture
+ f.write(self.ReprTable(**kargs) + '\n')
+ # fermeture
+ if kargs.get('FICHIER')<>None:
+ f.close()
+
+ def ReprTable(self,FORMAT='TABLEAU',dform=DicForm,**ignore):
+ """Représentation d'une Table ou d'une Colonne sous forme d'un tableau.
+ """
+ rows=self.rows
+ para=self.para
+ typ =self.type
+ if not type(para) in EnumTypes:
+ para=[self.para,]
+ typ =[self.type,]
+ # est-ce que l'attribut .type est renseigné ?
+ typdef=typ<>[None]*len(typ)
+ txt=[]
+ # ['']+ pour ajouter un séparateur en début de ligne
+ lspa=['',]
+ # lmax : largeur max des colonnes = max(form{K,R,I},len(parametre))
+ lmax=[]
+ for p in para:
+ t=typ[para.index(p)]
+ larg_max=max([len(str(p))] + \
+ [len(FMT(dform,k,t) % 0) for k in ('formK','formR','formI')])
+ lspa.append(FMT(dform,'formK',t,larg_max,str(p)) % p)
+ lmax.append(larg_max)
+ if typdef:
+ stype=dform['csep'].join([''] + \
+ [FMT(dform,'formK',typ[i],lmax[i]) % typ[i] for i in range(len(para))])
+ txt.append('')
+ txt.append('-'*80)
+ txt.append('')
+ ASTER=(FORMAT=='ASTER')
+ if ASTER:
+ txt.append('#DEBUT_TABLE')
+ if self.titr:
+ if ASTER:
+ txt.extend(['#TITRE '+lig for lig in self.titr.split('\n')])
+ else:
+ txt.append(self.titr)
+ txt.append(dform['csep'].join(lspa))
+ if ASTER and typdef:
+ txt.append(stype)
+ for r in rows:
+ lig=['']
+ empty=True
+ for v in para:
+ i=para.index(v)
+ t=typ[i]
+ rep=r.get(v,None)
+ if type(rep) is FloatType:
+ lig.append(FMT(dform,'formR',t,lmax[i]) % rep)
+ empty=False
+ elif type(rep) is IntType:
+ lig.append(FMT(dform,'formI',t,lmax[i]) % rep)
+ empty=False
+ else:
+ if rep==None:
+ rep='-'
+ else:
+ empty=False
+ s=FMT(dform,'formK',t,lmax[i],rep) % str(rep)
+ # format AGRAF = TABLEAU + '\' devant les chaines de caractères !
+ if FORMAT=='AGRAF':
+ s='\\'+s
+ lig.append(s)
+ if not empty:
+ txt.append(dform['csep'].join(lig))
+ if ASTER:
+ txt.append('#FIN_TABLE')
+ return dform['cfin'].join(txt)
+# ------------------------------------------------------------------------------
+ def ImprTabCroise(self,**kargs):
+ """Impression au format TABLEAU_CROISE d'une table ayant 3 paramètres.
+ """
+ # création du tableau croisé et impression au format TABLEAU
+ tabc=self.Croise()
+ kargs['FORMAT']='TABLEAU'
+ tabc.Impr(**kargs)
+# ------------------------------------------------------------------------------
+ def ImprGraph(self,**kargs):
+ """Impression au format XMGRACE : via le module Graph
+ """
+ args=kargs.copy()
+ if len(self.para)<>2:
+ UTMESS('A','Table','La table doit avoir exactement deux paramètres.')
+ return
+ lx, ly = [[v for v in getattr(self,p).values() if v<>None] for p in self.para]
+ # objet Graph
+ graph=Graph.Graph()
+ dicC={
+ 'Val' : [lx, ly],
+ 'Lab' : self.para,
+ }
+ if args['LEGENDE']==None: del args['LEGENDE']
+ Graph.AjoutParaCourbe(dicC, args)
+ graph.AjoutCourbe(**dicC)
+
+ # Surcharge des propriétés du graphique et des axes
+ # (bloc quasiment identique dans impr_fonction_ops)
+ if args.get('TITRE'): graph.Titre=args['TITRE']
+ if args.get('BORNE_X'):
+ graph.Min_X=args['BORNE_X'][0]
+ graph.Max_X=args['BORNE_X'][1]
+ if args.get('BORNE_Y'):
+ graph.Min_Y=args['BORNE_Y'][0]
+ graph.Max_Y=args['BORNE_Y'][1]
+ if args.get('LEGENDE_X'): graph.Legende_X=args['LEGENDE_X']
+ if args.get('LEGENDE_Y'): graph.Legende_Y=args['LEGENDE_Y']
+ if args.get('ECHELLE_X'): graph.Echelle_X=args['ECHELLE_X']
+ if args.get('ECHELLE_Y'): graph.Echelle_Y=args['ECHELLE_Y']
+ if args.get('GRILLE_X'): graph.Grille_X=args['GRILLE_X']
+ if args.get('GRILLE_Y'): graph.Grille_Y=args['GRILLE_Y']
+
+ try:
+ graph.Trace(**args)
+ except TypeError:
+ UTMESS('A','Table','Les cellules ne doivent contenir que des nombres réels')
+
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+class Table(TableBase):
+ """Une table est construite comme une liste de lignes, chaque ligne est
+ un dictionnaire.
+ On crée puis on ajoute les lignes avec la méthode append :
+ t=Table()
+ t.append(dict(a=1,b=2))
+ t.append(dict(a=3,b=4))
+ La méthode __iter__ définit un itérateur sur les lignes de la table,
+ __repr__ retourne une représentation de la table, utilisée par "print t".
+ Grace à la classe Colonne et à sa méthode _extract, il est possible
+ de construire une sous-table qui satisfait un critère donné.
+ Le critère est donné par une fonction Python qui retourne vrai
+ ou faux si la valeur d'une colonne respecte le critère ou non.
+ Exemple:
+ def critere(valeur):
+ return valeur < 10
+ soustable = t.a._extract(critere)
+ t.a retourne un objet intermédiaire de la classe Colonne qui mémorise
+ le nom de la colonne demandée (a, ici).
+ """
+ def __init__(self, rows=[], para=[], typ=[], titr=''):
+ """Constructeur de la Table :
+ rows : liste des lignes (dict)
+ para : liste des paramètres
+ type : liste des types des paramètres
+ titr : titre de la table
+ """
+ self.rows=[r for r in rows if r.values()<>[None]*len(r.values())]
+ self.para=list(para)
+ if len(typ)==len(self.para):
+ self.type=list(typ)
+ else:
+ self.type=[None]*len(self.para)
+ self.titr=titr
+
+ def append(self, obj):
+ """Ajoute une ligne (type dict) à la Table"""
+ self.rows.append(obj)
+
+ def __iter__(self):
+ """Itère sur les lignes de la Table"""
+ return iter(self.rows)
+
+ def __getattr__(self, column):
+ """Construit un objet intermediaire (couple table, colonne)"""
+ typ=None
+ if not column in self.para:
+ column=''
+ else:
+ typ=self.type[self.para.index(column)]
+ return Colonne(self, column, typ)
+
+ def sort(self, CLES=None, ORDRE='CROISSANT'):
+ """Tri de la table.
+ CLES : liste des clés de tri
+ ORDRE : CROISSANT ou DECROISSANT (de longueur 1 ou len(keys))
+ """
+ # par défaut, on prend tous les paramètres
+ if CLES==None:
+ CLES=self.para[:]
+ if not type(CLES) in EnumTypes:
+ CLES=[CLES,]
+ else:
+ CLES=list(CLES)
+ self.rows=sort_table(self.rows, self.para, CLES, (ORDRE=='DECROISSANT'))
+# if not type(order) in EnumTypes:
+# order=[order,]
+# print 'TRI clés=%s, order=%s' % (keys,order)
+# # on ne garde que le premier si les longueurs sont différentes
+# if len(order)<>len(keys):
+# order=[order[0],]
+# else:
+# # si toutes les valeurs sont identiques, on peut ne garder que la 1ère
+# d={}
+# for o in order: d[o]=None
+# if len(order)<>len(keys) or len(d.keys())==1:
+# order=[order[0],]
+# if len(order)==1:
+# self.rows=sort_table(self.rows, self.para, keys, (order[0]=='DECROISSANT'))
+# else:
+# # de la dernière clé à la première
+# for k,o in [(keys[i],order[i]) for i in range(len(keys)-1,-1,-1)]:
+# print 'TRI : clé=%s, order=%s' % (k,o)
+# self.rows=sort_table(self.rows, self.para, [k], (o=='DECROISSANT'))
+
+ def __delitem__(self, args):
+ """Supprime les colonnes correspondantes aux éléments de args """
+ if not type(args) in EnumTypes:
+ args=[args,]
+ new_rows=self.rows
+ new_para=self.para
+ new_type=self.type
+ for item in args:
+ del new_type[new_para.index(item)]
+ new_para.remove(item)
+ for line in new_rows : del line[item]
+ return Table(new_rows, new_para, new_type, self.titr)
+
+ def __getitem__(self, args):
+ """Extrait la sous table composée des colonnes dont les paramètres sont dans args """
+ if not type(args) in EnumTypes:
+ args=[args,]
+ else:
+ args=list(args)
+ #print '<getitem> args=',args
+ new_rows=[]
+ new_para=args
+ new_type=[]
+ for item in new_para:
+ if not item in self.para:
+ return Table()
+ new_type.append(self.type[self.para.index(item)])
+ for line in self:
+ new_line={}
+ for item in new_para:
+ new_line[item]=line.get(item)
+ new_rows.append(new_line)
+ return Table(new_rows, new_para, new_type, self.titr)
+
+ def __and__(self, other):
+ """Intersection de deux tables (opérateur &)"""
+ if other.para<>self.para:
+ UTMESS('A','Table','Les paramètres sont différents')
+ return Table()
+ else:
+ tmp = [ r for r in self if r in other.rows ]
+ return Table(tmp, self.para, self.type, self.titr)
+
+ def __or__(self, other):
+ """Union de deux tables (opérateur |)"""
+ if other.para<>self.para:
+ UTMESS('A','Table','Les paramètres sont différents')
+ return Table()
+ else:
+ tmp = self.rows[:]
+ tmp.extend([ r for r in other if r not in self ])
+ return Table(tmp, self.para, self.type[:], self.titr)
+
+ def values(self):
+ """Renvoie la table sous la forme d'un dictionnaire de listes dont les
+ clés sont les paramètres.
+ """
+ dico={}
+ for column in self.para:
+ dico[column]=Colonne(self, column).values()
+ return dico
+
+ def Array(self,Para,Champ):
+ """Renvoie sous forme de NumArray le résultat d'une extraction dans une table
+ méthode utile à macr_recal
+ """
+ import Numeric
+ __Rep = self[Para,Champ].values()
+ F=Numeric.zeros((len(__Rep[Para]),2),Numeric.Float)
+ for i in range(len(__Rep[Para])):
+ F[i][0] = __Rep[Para][i]
+ F[i][1] = __Rep[Champ][i]
+ del(__Rep)
+ return F
+
+ def Croise(self):
+ """Retourne un tableau croisé P3(P1,P2) à partir d'une table ayant
+ trois paramètres (P1, P2, P3).
+ """
+ if len(self.para)<>3:
+ UTMESS('A','Table','La table doit avoir exactement trois paramètres.')
+ return Table()
+ py, px, pz = self.para
+ ly, lx, lz = [getattr(self,p).values() for p in self.para]
+ new_rows=[]
+ #lpz='%s=f(%s,%s)' % (pz,px,py)
+ lpz='%s/%s' % (px,py)
+ new_para=[lpz,]
+ # attention aux doublons dans lx et ly
+ for it in ly:
+ if it<>None and new_para.count(it)==0:
+ new_para.append(it)
+ newx=[]
+ for it in lx:
+ if it<>None and newx.count(it)==0:
+ newx.append(it)
+ for x in newx:
+ if x<>None:
+ d={ lpz : x, }
+ taux = (getattr(self,px)==x)
+ for dz in taux.rows:
+ d[dz[py]]=dz[pz]
+ new_rows.append(d)
+ new_type=[self.type[0],] + [self.type[2]]*len(ly)
+ new_titr=self.titr
+ if new_titr<>'': new_titr+='\n'
+ new_titr+=pz + ' FONCTION DE ' + px + ' ET ' + py
+ return Table(new_rows, new_para, new_type, new_titr)
+
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+class Colonne(TableBase):
+ """Classe intermédiaire pour mémoriser un couple (table, nom de colonne)
+ et exprimer les critères d'extraction sous une forme naturelle en python
+ en surchargeant les operateurs <, >, <> et =.
+ Alors on peut écrire la requete simple :
+ soustable=t.a<10
+ Ainsi que des requetes plus complexes :
+ soustable=t.a<10 & t.b <4
+ ou
+ soustable=t.a<10 | t.b <4
+ Les "alias" EQ, NE, LE, LT, GE, GT permettent à la macro IMPR_TABLE
+ d'utiliser directement le mot-clé utilisateur CRIT_COMP défini dans le
+ catalogue : getattr(Table,CRIT_COMP).
+ """
+ def __init__(self, table, column, typ=None):
+ """Constructeur (objet Table associé, paramètre de la colonne, type du
+ paramètre).
+ """
+ self.Table=table
+ self.rows=self.Table.rows
+ self.para=column
+ self.type=typ
+ self.titr=''
+
+ def _extract(self, fun):
+ """Construit une table avec les lignes de self.Table
+ dont l'élément de nom self.para satisfait le critère fun,
+ fun est une fonction qui retourne vrai ou faux
+ """
+ return Table([row for row in self.Table if fun(row.get(self.para))], self.Table.para, self.Table.type, self.Table.titr)
+
+ def __le__(self, VALE):
+ return self._extract(lambda v: v<>None and v<=VALE)
+
+ def __lt__(self, VALE):
+ return self._extract(lambda v: v<>None and v<VALE)
+
+ def __ge__(self, VALE):
+ return self._extract(lambda v: v<>None and v>=VALE)
+
+ def __gt__(self, VALE):
+ return self._extract(lambda v: v<>None and v>VALE)
+
+ def __eq__(self, VALE, CRITERE='RELATIF', PRECISION=0.):
+ if type(VALE) in EnumTypes :
+ return self._extract(lambda v: v in VALE)
+ if PRECISION==0. or not type(VALE) in NumberTypes:
+ if type(VALE) in StringTypes:
+ return self._extract(lambda v: v<>None and str(v).strip()==VALE.strip())
+ else:
+ return self._extract(lambda v: v==VALE)
+ else:
+ if CRITERE=='ABSOLU':
+ vmin=VALE-PRECISION
+ vmax=VALE+PRECISION
+ else:
+ vmin=(1.-PRECISION)*VALE
+ vmax=(1.+PRECISION)*VALE
+ return self._extract(lambda v: v<>None and vmin<v<vmax)
+
+ def __ne__(self, VALE, CRITERE='RELATIF', PRECISION=0.):
+ if type(VALE) in EnumTypes :
+ return self._extract(lambda v: v not in VALE)
+ if PRECISION==0. or not type(VALE) in NumberTypes:
+ if type(VALE) in StringTypes:
+ return self._extract(lambda v: v<>None and str(v).strip()<>VALE.strip())
+ else:
+ return self._extract(lambda v: v<>VALE)
+ else:
+ if CRITERE=='ABSOLU':
+ vmin=VALE-PRECISION
+ vmax=VALE+PRECISION
+ else:
+ vmin=(1.-PRECISION)*VALE
+ vmax=(1.+PRECISION)*VALE
+ return self._extract(lambda v: v<>None and (v<vmin or vmax<v))
+
+ def MAXI(self):
+ # important pour les performances de récupérer le max une fois pour toutes
+ maxi=max(self)
+ return self._extract(lambda v: v==maxi)
+
+ def MINI(self):
+ # important pour les performances de récupérer le min une fois pour toutes
+ mini=min(self)
+ return self._extract(lambda v: v==mini)
+
+ def ABS_MAXI(self):
+ # important pour les performances de récupérer le max une fois pour toutes
+ abs_maxi=max([abs(v) for v in self.values() if type(v) in NumberTypes])
+ return self._extract(lambda v: v==abs_maxi or v==-abs_maxi)
+
+ def ABS_MINI(self):
+ # important pour les performances de récupérer le min une fois pour toutes
+ abs_mini=min([abs(v) for v in self.values() if type(v) in NumberTypes])
+ # tester le type de v est trop long donc pas de abs(v)
+ return self._extract(lambda v: v==abs_mini or v==-abs_mini)
+
+ def __iter__(self):
+ """Itère sur les éléments de la colonne"""
+ for row in self.Table:
+ # si l'élément n'est pas présent on retourne None
+ yield row.get(self.para)
+ #yield row[self.para]
+
+ def __getitem__(self, i):
+ """Retourne la ième valeur d'une colonne"""
+ return self.values()[i]
+
+ def values(self):
+ """Renvoie la liste des valeurs"""
+ return [r[self.para] for r in self.Table]
+
+ # équivalences avec les opérateurs dans Aster
+ LE=__le__
+ LT=__lt__
+ GE=__ge__
+ GT=__gt__
+ EQ=__eq__
+ NE=__ne__
+ def VIDE(self) : return self.__eq__(None)
+ def NON_VIDE(self): return self.__ne__(None)
+
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+def sort_table(rows,l_para,w_para,reverse=False):
+ """Sort list of dict.
+ rows : list of dict
+ l_para : list of the keys of dict
+ w_para : keys of the sort
+ """
+ c_para=[i for i in l_para if i not in w_para]
+ new_rows=rows
+ for i in w_para :
+ new_key= '__'+str(w_para.index(i))+i
+ for row in new_rows :
+ row[new_key]=row[i]
+ del row[i]
+ for i in c_para :
+ new_key= '___'+i
+ for row in new_rows :
+ row[new_key]=row[i]
+ del row[i]
+ new_rows.sort()
+ if reverse:
+ new_rows.reverse()
+ for i in w_para :
+ old_key= '__'+str(w_para.index(i))+i
+ for row in new_rows :
+ row[i]=row[old_key]
+ del row[old_key]
+ for i in c_para :
+ old_key= '___'+i
+ for row in new_rows :
+ row[i]=row[old_key]
+ del row[old_key]
+ return new_rows
+
+# ------------------------------------------------------------------------------
+def FMT(dform, nform, typAster=None, larg=0, val=''):
+ """Retourne un format d'impression Python à partir d'un type Aster ('R','I',
+ 'K8', 'K16'...). Si typAster==None, retourne dform[nform].
+ larg : largeur minimale du format (val permet de ne pas ajouter des blancs
+ si la chaine à afficher est plus longue que le format, on prend le partie
+ de ne pas tronquer les chaines)
+ """
+ if typAster==None:
+ fmt=dform[nform]
+ elif typAster in ('I', 'R'):
+ if nform=='formK':
+ # convertit %12.5E en %-12s
+ fmt=re.sub('([0-9]+)[\.0-9]*[diueEfFgG]+','-\g<1>s',dform['form'+typAster])
+ #print nform, typAster, fmt
+ else:
+ fmt=dform[nform]
+ else:
+ # typAster = Kn
+ fmt='%-'+typAster[1:]+'s'
+ # on ajoute éventuellement des blancs pour atteindre la largeur demandée
+ if larg<>0:
+ fmt=' '*max(min(larg-len(val),larg-len(fmt % 0)),0) + fmt
+ return fmt
+
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+# ------------------------------------------------------------------------------
+if __name__ == "__main__":
+ listdic = [
+ {'NOEUD': 'N1' ,'NUME_ORDRE': 1 ,'INST': 0.5, 'DX': -0.00233, 'COOR_Y': 0.53033,},
+ {'NOEUD': 'N1' ,'NUME_ORDRE': 2 ,'INST': 1.0, 'DX': -0.00467, 'COOR_Y': 0.53033,},
+ {'NOEUD': 'N1' ,'NUME_ORDRE': 3 ,'INST': 1.5, 'DX': -0.00701, 'COOR_Y': 0.53033,},
+ {'NOEUD': 'N1' ,'NUME_ORDRE': 4 ,'INST': 2.0, 'DX': -0.00934, 'COOR_Y': 0.53033,},
+ {'NOEUD': 'N1' ,'NUME_ORDRE': 5 ,'INST': 2.5, 'DX': -0.01168, 'COOR_Y': 0.53033,},
+ {'NOEUD': 'N2' ,'NUME_ORDRE': 11,'INST': 5.5, 'DX': -0.00233, 'COOR_Y': 0.53033,},
+ {'NOEUD': 'N2' ,'NUME_ORDRE': 12,'INST': 6.0, 'DX': -0.00467, 'COOR_Y': 0.53033,},
+ {'NOEUD': 'N2' ,'NUME_ORDRE': 13,'INST': 6.5, 'DX': -0.00701, 'COOR_Y': 0.53033,},
+ {'NOEUD': 'N2' ,'NUME_ORDRE': 14,'INST': 7.0, 'DX': -0.00934, 'COOR_Y': 0.53033,},
+ {'NOEUD': 'N2' ,'NUME_ORDRE': 15,'INST': 7.5, 'DX': -0.01168, 'COOR_Y': 0.53033,},
+ ]
+ import random
+ random.shuffle(listdic)
+ listpara=['NOEUD','NUME_ORDRE','INST','COOR_Y','DX']
+ listtype=['K8','I','R','R','R']
+ t=Table(listdic,listpara,listtype)
+
+ tb=t[('NOEUD','DX')]
+ print tb.para
+ print tb.type
+
+ print
+ print "------Table initiale----"
+ print t
+ print
+ print "--------- CRIT --------"
+ print t.NUME_ORDRE <=5
+ print
+ print "------- CRIT & CRIT -----"
+ print (t.NUME_ORDRE < 10) & (t.INST >=1.5)
+ print
+ print "----- EQ maxi / min(col), max(col) ------"
+ print t.DX == max(t.DX)
+ print min(t.DX)
+ print max(t.DX)
+ print "------ getitem sur 2 paramètres ------"
+ print t.NUME_ORDRE
+ print t.DX
+ print t['DX','NUME_ORDRE']
+ print "------ sort sur INST ------"
+ t.sort('INST')
+ print t
+
+ print "------- TABLEAU_CROISE ------"
+ tabc=t['NOEUD','INST','DX']
+ tabc.Impr(FORMAT='TABLEAU_CROISE')
+
+ N=5
+ ldic=[]
+ for i in range(N):
+ ldic.append({'IND':float(i), 'VAL' : random.random()*i})
+ para=['IND','VAL']
+ t3=Table(ldic, para, titr='Table aléatoire')
+ col=t3.VAL.ABS_MAXI()
+ col=t3.VAL.MINI()
+
+ t3.sort('VAL','IND')
+
+ tg=tabc['INST','DX'].DX.NON_VIDE()
+ #tg.Impr(FORMAT='XMGRACE')
+
+ g=Graph.Graph()
+ g.Titre="Tracé d'une fonction au format TABLEAU"
+ g.AjoutCourbe(Val=[tg.INST.values(), tg.DX.values()], Lab=['INST','DX'])
+ g.Trace(FORMAT='TABLEAU')
+
+# t.Impr(PAGINATION='NOEUD')
+ t.Impr(PAGINATION=('NOEUD','INST'))
+
--- /dev/null
+#@ MODIF UniteAster Utilitai DATE 11/05/2005 AUTEUR MCOURTOI M.COURTOIS
+# -*- coding: iso-8859-1 -*-
+# CONFIGURATION MANAGEMENT OF EDF VERSION
+# ======================================================================
+# COPYRIGHT (C) 1991 - 2005 EDF R&D WWW.CODE-ASTER.ORG
+# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
+# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
+# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
+# (AT YOUR OPTION) ANY LATER VERSION.
+#
+# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
+# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
+# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
+# GENERAL PUBLIC LICENSE FOR MORE DETAILS.
+#
+# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
+# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
+# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
+# ======================================================================
+
+import types
+
+import aster
+from Cata.cata import _F, DEFI_FICHIER, INFO_EXEC_ASTER, DETRUIRE
+
+#-------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------
+class UniteAster:
+ """Classe pour manipuler les fichiers en Python en accord avec les unités
+ logiques utilisées en Fortran.
+ De manière analogue au Fortran, les états possibles d'une unité sont :
+ 'F' : fermé, 'O' : ouvert, 'R' : réservé.
+
+ Méthodes :
+ Nom : Retourne le nom du fichier associé à une unité,
+ Etat : Retourne l'état d'une unité,
+ Libre : Retourne un numéro d'unité libre,
+ EtatInit : Remet une, plusieurs ou toutes les unités dans leur état initial.
+
+ Méthode privée :
+ _setinfo : pour remplir le dictionnaire des 'infos'
+ Attribut privé :
+ infos[numéro unité] = { 'nom' : x, 'etat' : x , 'etat_init' : x }
+ """
+#-------------------------------------------------------------------------------
+ def __init__(self):
+ """Initialise le dictionnaire des unités.
+ """
+ self.infos = {}
+
+#-------------------------------------------------------------------------------
+ def _setinfo(self, ul):
+ """Remplit les infos de l'unité 'ul'.
+ """
+ # ul peut etre un entier Aster
+ try:
+ unit = ul.valeur
+ except:
+ unit = int(ul)
+ # Si la clé n'existe pas
+ ini = False
+ if not self.infos.has_key(unit):
+ self.infos[unit] = {}
+ self.infos[unit]['nom'] = ''
+ self.infos[unit]['etat'] = '?'
+ self.infos[unit]['etat_init'] = '?'
+ ini = True
+
+ __tab=INFO_EXEC_ASTER(UNITE=unit, LISTE_INFO=('ETAT_UNITE'))
+
+ # O:ouvert, F:fermé, R:réservé
+ self.infos[unit]['etat'] = __tab['ETAT_UNITE',1].strip()[0]
+ if ini:
+ self.infos[unit]['etat_init'] = self.infos[unit]['etat']
+
+ # nom du fichier
+ if self.infos[unit]['etat'] in ['O', 'R']:
+ nomfich=''.join([__tab['NOMFIC%d' % i,1] for i in range(1,5)]).strip()
+ elif self.infos[unit]['etat'] == 'F':
+ nomfich='fort.'+str(unit)
+ else:
+ message = "Etat de l'unité inconnu : %s" % self.infos[unit]['etat']
+ print __tab.EXTR_TABLE()
+ raise aster.FatalError,"<F> <UniteAster._setinfo> %s" % message
+ self.infos[unit]['nom'] = nomfich
+ DETRUIRE(CONCEPT=_F(NOM=__tab))
+
+#-------------------------------------------------------------------------------
+ def Libre(self, nom=None):
+ """Réserve et retourne une unité libre en y associant, s'il est fourni,
+ le fichier 'nom'.
+ """
+ __tab=INFO_EXEC_ASTER(LISTE_INFO=('UNITE_LIBRE'))
+ unit = __tab['UNITE_LIBRE',1]
+ DETRUIRE(CONCEPT=_F(NOM=__tab))
+ if nom==None:
+ nom='fort.'+str(unit)
+
+ # Si la clé existe, c'est que le fichier n'était pas libre
+ if self.infos.has_key(unit):
+ message = "Cette unité est déjà affectée au fichier %s" % \
+ self.infos[unit]['nom']
+ raise aster.FatalError,"<F> <UniteAster.Libre> %s" % message
+
+ DEFI_FICHIER(ACTION='RESERVER', UNITE=unit , FICHIER=nom.strip())
+ self.infos[unit] = {}
+ self.infos[unit]['nom'] = nom.strip()
+ self.infos[unit]['etat'] = 'R'
+ self.infos[unit]['etat_init'] = 'F'
+ return unit
+
+#-------------------------------------------------------------------------------
+ def Nom(self, ul):
+ """Retourne le nom du fichier associé à l'unité 'ul'.
+ """
+ # ul peut etre un entier Aster
+ try:
+ unit = ul.valeur
+ except:
+ unit = int(ul)
+ # Si la clé n'existe pas
+ if not self.infos.has_key(unit):
+ self._setinfo(unit)
+ return self.infos[unit]['nom']
+
+#-------------------------------------------------------------------------------
+ def Etat(self, ul, **kargs):
+ """Retourne l'état de l'unité si 'etat' n'est pas fourni
+ et/ou change son état :
+ kargs['etat'] : nouvel état,
+ kargs['TYPE'] : type du fichier à ouvrir ASCII/BINARY/LIBRE,
+ kargs['ACCES'] : type d'accès NEW/APPEND/OLD.
+ """
+ # ul peut etre un entier Aster
+ try:
+ unit = ul.valeur
+ except:
+ unit = int(ul)
+ # Si la clé n'existe pas
+ if not self.infos.has_key(unit):
+ self._setinfo(unit)
+ if not kargs.has_key('etat'):
+ return self.infos[unit]['etat']
+
+ # En fonction de la demande, on bascule son état ou pas
+ new = kargs.get('etat')
+ if not new in ['R', 'F', 'O']:
+ message = "Nouvel état de l'unité incorrect : %s" % new
+ raise aster.FatalError,"<F> <UniteAster.Etat> %s" % message
+
+ if self.infos[unit]['etat'] == new:
+ pass
+ elif new == 'R':
+ if self.infos[unit]['etat'] == 'O':
+ DEFI_FICHIER(ACTION='LIBERER', UNITE=unit)
+ DEFI_FICHIER(ACTION = 'RESERVER',
+ UNITE = unit,
+ FICHIER = self.infos[unit]['nom'])
+ elif new == 'F':
+ DEFI_FICHIER(ACTION='LIBERER', UNITE=unit)
+ elif new == 'O':
+ if self.infos[unit]['etat'] == 'R':
+ DEFI_FICHIER(ACTION='LIBERER', UNITE=unit)
+ DEFI_FICHIER(ACTION ='ASSOCIER',
+ UNITE = unit,
+ FICHIER = self.infos[unit]['nom'],
+ TYPE = kargs.get('TYPE', 'ASCII'),
+ ACCES = kargs.get('ACCES', 'APPEND'),)
+ self.infos[unit]['etat'] = new
+ return self.infos[unit]['etat']
+
+#-------------------------------------------------------------------------------
+ def EtatInit(self, ul=None):
+ """Remet l'unité 'ul' dans son état initial.
+ Si 'ul' est omis, toutes les unités sont remises dans leur état initial.
+ """
+ if ul == None:
+ for uli, vul in self.infos.items():
+ self.Etat(uli, etat=vul['etat_init'])
+ else:
+ if not type(ul) in [types.ListType, types.TupleType]:
+ ul=[ul,]
+ for u in ul:
+ # u peut etre un entier Aster
+ try:
+ unit = u.valeur
+ except:
+ unit = int(u)
+ # Si la clé n'existe pas
+ if not self.infos.has_key(unit):
+ self._setinfo(unit)
+ else:
+ self.Etat(unit, etat=self.infos[unit]['etat_init'])
--- /dev/null
+#@ MODIF Utmess Utilitai DATE 30/11/2004 AUTEUR MCOURTOI M.COURTOIS
+# -*- coding: iso-8859-1 -*-
+# CONFIGURATION MANAGEMENT OF EDF VERSION
+# ======================================================================
+# COPYRIGHT (C) 1991 - 2004 EDF R&D WWW.CODE-ASTER.ORG
+# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
+# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
+# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
+# (AT YOUR OPTION) ANY LATER VERSION.
+#
+# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
+# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
+# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
+# GENERAL PUBLIC LICENSE FOR MORE DETAILS.
+#
+# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
+# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
+# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
+# ======================================================================
+
+import sys
+import aster
+
+def UTMESS(code, sprg, texte):
+ """Utilitaire analogue à la routine fortran UTMESS.
+ code : 'A', 'E', 'S', 'F'
+ sprg : nom du module, classe ou fonction python où l'on se trouve
+ texte : contenu du message
+ """
+ fmt='\n <%s> <%s> %s\n\n'
+ UL={
+ 'MESSAGE' : 6,
+ 'RESULTAT' : 8,
+ #'ERREUR' : 9,
+ }
+ # On importe la définition des commandes à utiliser dans la macro
+# if jdc:
+# DEFI_FICHIER = jdc.get_cmd('DEFI_FICHIER')
+# else:
+# # on se limite au print !
+# UL={ 'MESSAGE' : 6, }
+ try:
+ from Cata.cata import DEFI_FICHIER
+ except ImportError:
+ # on se limite au print !
+ UL={ 'MESSAGE' : 6, }
+
+ reason=fmt % (code, sprg, texte)
+
+ for nom,ul in UL.items():
+ if ul<>6:
+ DEFI_FICHIER(ACTION='LIBERER', UNITE=ul, )
+ f=open('fort.'+str(ul),'a')
+ else:
+ f=sys.stdout
+ # écriture du message
+ f.write(reason)
+
+ if ul<>6:
+ f.close()
+ DEFI_FICHIER(ACTION='ASSOCIER', UNITE=ul, TYPE='ASCII', ACCES='APPEND')
+
+ if code=='S':
+ raise aster.error, reason
+ elif code=='F':
+ raise aster.FatalError, reason
--- /dev/null
+#@ MODIF __init__ Utilitai DATE 20/09/2004 AUTEUR DURAND C.DURAND
+# -*- coding: iso-8859-1 -*-
+# CONFIGURATION MANAGEMENT OF EDF VERSION
+# ======================================================================
+# COPYRIGHT (C) 1991 - 2003 EDF R&D WWW.CODE-ASTER.ORG
+# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
+# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
+# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
+# (AT YOUR OPTION) ANY LATER VERSION.
+#
+# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
+# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
+# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
+# GENERAL PUBLIC LICENSE FOR MORE DETAILS.
+#
+# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
+# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
+# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
+# ======================================================================
+
+
--- /dev/null
+#@ MODIF courbes Utilitai DATE 14/09/2004 AUTEUR MCOURTOI M.COURTOIS
+# -*- coding: iso-8859-1 -*-
+# CONFIGURATION MANAGEMENT OF EDF VERSION
+# ======================================================================
+# COPYRIGHT (C) 1991 - 2004 EDF R&D WWW.CODE-ASTER.ORG
+# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
+# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
+# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
+# (AT YOUR OPTION) ANY LATER VERSION.
+#
+# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
+# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
+# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
+# GENERAL PUBLIC LICENSE FOR MORE DETAILS.
+#
+# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
+# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
+# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
+# ======================================================================
+
+#==================================================
+# fonction "COURBES"
+# usage : permet de tracer des courbes en interactif
+# avec XMGRACE ou dans un fichier postscript
+#==================================================
+
+import Stanley
+from Stanley import xmgrace
+from Stanley import as_courbes
+
+def COURBES(listcourb,titre=' ',soustitre=' ',legx=' ',legy=' ',bornex=None,borney=None,fichier=None):
+
+# ARGUMENTS
+
+# listcourb : tuple de courbes, chaque courbe etant definie soit par
+# (TABLE1, NOM_PARA_X, TABLE2, NOM_PARA_Y, LEGENDE)
+# soit par :
+# (FONCTION,LEGENDE)
+# titre et sous_titre : facultatifs, titre et sous-tritre du graphique
+# legx, legy : facultatifs, legendes des axes
+# bornex, borney : facultatifs, bornes sur les axes
+# fichier : facultatif : sortie au format postscript si present
+#
+# exemples d'appel :
+#--------------------
+# courb1=(SYYPRT,'ABSC_CURV',SYYPRT,'SIYY','PRT')
+# courb2=(SYYMLC10,'ABSC_CURV',SYYMLC10,'SIYY','MLC10')
+# courb3=(SYYML100,'ABSC_CURV',SYYML100,'SIYY','MLC100')
+# listcourb=(courb1,courb2,courb3)
+# COURBES(listcourb,titre='Plaque trouee',legx='Abcisses curvilignes',legy='Contraintes (MPa)',bornex=(0,100),borney=(500,1000))
+# fonc1=(F1,'F_PRT')
+# fonc2=(F2,'F_MLC10')
+# fonc3=(F3,'F_MLC100')
+# listfonc=(fonc1,fonc2,fonc3)
+# COURBES(listfonc,titre='Fonctions')
+# postscript
+# COURBES(listfonc,titre='Plaque trouee',fichier='./fort.24')
+#--------------------------------------------------------------
+
+# initialisation du trace de courbes
+
+ if (fichier!=None):
+ graphe=xmgrace.Xmgr(10,' -hardcopy -nosafe')
+ print "Nombre de courbes ",len(listcourb)," sur le fichier :",fichier
+
+ else:
+ graphe=xmgrace.Xmgr(10,' -noask')
+ print "Nombre de courbes ",len(listcourb)
+
+ graphe.Nouveau_graphe()
+
+# dimensionnement des axes
+ if bornex != None :
+ xmin=list(bornex)[0]
+ xmax=list(bornex)[1]
+ ctest1 = as_courbes.Courbe()
+ ctest1.x=[xmin,xmax]
+ ctest1.y=[0.0,0.0]
+ graphe.Courbe(ctest1)
+
+ if borney != None :
+ ymin=list(borney)[0]
+ ymax=list(borney)[1]
+ ctest2 = as_courbes.Courbe()
+ ctest2.x=[0.0,0.0]
+ ctest2.y=[ymin,ymax]
+ graphe.Courbe(ctest2)
+
+ if titre != None :
+ if soustitre != None :
+ graphe.Titre(titre,soustitre)
+ else :
+ graphe.Titre(titre,' ')
+
+ if legx != None :
+ graphe.Axe_x(legx)
+
+ if legy != None :
+ graphe.Axe_y(legy)
+
+ k = 0
+
+ for courbi in listcourb:
+ sigi = as_courbes.Courbe()
+
+ try :
+ # cas d une table
+ sigi.Lire_x(courbi[0],courbi[1])
+ sigi.Lire_y(courbi[2],courbi[3])
+ legende=courbi[4]
+ except :
+ # cas d une fonction
+ sigi.x,sigi.y=courbi[0].Valeurs()
+ legende=courbi[1]
+
+ graphe.Courbe(sigi,legende)
+ graphe.Send('WITH G'+repr(graphe.gr_act))
+ graphe.Send('S' + str(k) + ' SYMBOL ' + str(k+2))
+ graphe.Send('S' + str(k) + ' SYMBOL SIZE 0.5')
+ graphe.Send('S' + str(k) + ' SYMBOL COLOR '+str(k+2))
+ graphe.Send('S' + str(k) + ' LINE COLOR '+str(k+2))
+ k = k + 1
+ graphe.Send('REDRAW')
+
+ if (fichier!=None):
+ graphe.Sortie_EPS(fichier)
+ graphe.Fermer()
+ else:
+ graphe.Attendre()
+
+ k=0
+
+#===========================================
+
+
--- /dev/null
+#@ MODIF funct_root Utilitai DATE 14/09/2004 AUTEUR MCOURTOI M.COURTOIS
+# -*- coding: iso-8859-1 -*-
+################################################################################
+# Mathematical utility routines
+# Copyright (C) 1999, Wesley Phoa
+#
+# Reference: Numerical Recipes in C
+# [[[[extraits]]]]
+
+class BracketingException(Exception):
+ pass
+
+class RootFindingException(Exception):
+ pass
+
+class MinimizationException(Exception):
+ pass
+
+GOLDEN = (1+5**.5)/2
+
+#
+# MISCELLANEOUS
+#
+
+def sgn(x):
+ if x==0:
+ return 0
+ else:
+ return x/abs(x)
+
+#
+# UNIVARIATE ROOT FINDING
+#
+
+def bracket_root(f, interval, max_iterations=50):
+ """\
+Given a univariate function f and a tuple interval=(x1,x2),
+return a new tuple (bracket, fnvals) where bracket=(x1,x2)
+brackets a root of f and fnvals=(f(x1),f(x2)).
+ """
+ x1, x2 = interval
+ if x1==x2:
+ raise BracketingException("initial interval has zero width")
+ elif x2<x1:
+ x1, x2 = x2, x1
+ f1, f2 = f(x1), f(x2)
+ for j in range(max_iterations):
+ while f1*f2 >= 0: # not currently bracketed
+ if abs(f1)<abs(f2):
+ x1 = x1 + GOLDEN*(x1-x2)
+ else:
+ x2 = x2 + GOLDEN*(x2-x1)
+ f1, f2 = f(x1), f(x2)
+ return (x1, x2), (f1, f2)
+ raise BracketingException("too many iterations")
+
+def ridder_root(f, bracket, fnvals=None, accuracy=1e-6, max_iterations=50):
+ """\
+Given a univariate function f and a tuple bracket=(x1,x2) bracketing a root,
+find a root x of f using Ridder s method. Parameter fnvals=(f(x1),f(x2)) is optional.
+ """
+ x1, x2 = bracket
+ if fnvals==None:
+ f1, f2 = f(x1), f(x2)
+ else:
+ f1, f2 = fnvals
+ if f1==0:
+ return x1
+ elif f2==0:
+ return x2
+ elif f1*f2>=0:
+ raise BracketingException("initial interval does not bracket a root")
+ x4 = 123456789.
+ for j in range(max_iterations):
+ x3 = (x1+x2)/2
+ f3 = f(x3)
+ temp = f3*f3 - f1*f2
+ x4, x4old = x3 + (x3-x1)*sgn(f1-f2)*f3/temp**.5, x4
+ f4 = f(x4)
+ if f1*f4<0: # x1 and x4 bracket root
+ x2, f2 = x4, f4
+ else: # x4 and x2 bracket root
+ x1, f1 = x4, f4
+ if min(abs(x1-x2),abs(x4-x4old))<accuracy or temp==0:
+ return x4
+ raise RootFindingException("too many iterations")
+
+def root(f, interval=(0.,1.), accuracy=1e-4, max_iterations=50):
+ """\
+Given a univariate function f and an optional interval (x1,x2),
+find a root of f using bracket_root and ridder_root.
+ """
+ bracket, fnvals = bracket_root(f, interval, max_iterations)
+ return ridder_root(f, bracket, fnvals, accuracy, max_iterations)
+
--- /dev/null
+#@ MODIF partition Utilitai DATE 10/02/2005 AUTEUR ASSIRE A.ASSIRE
+# -*- coding: iso-8859-1 -*-
+# CONFIGURATION MANAGEMENT OF EDF VERSION
+# ======================================================================
+# COPYRIGHT (C) 1991 - 2004 EDF R&D WWW.CODE-ASTER.ORG
+# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
+# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
+# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
+# (AT YOUR OPTION) ANY LATER VERSION.
+#
+# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
+# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
+# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
+# GENERAL PUBLIC LICENSE FOR MORE DETAILS.
+#
+# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
+# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
+# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
+# ======================================================================
+# RESPONSABLE ASSIRE A.ASSIRE
+
+import aster
+import string, os, time, sys, UserList, types
+import Numeric as N
+
+from Accas import _F
+from Noyau.N_utils import AsType
+
+# ------------------------------------------------------------------ #
+
+def enleve_doublons_liste(liste):
+ """ A partir d une liste qui peut contenir des doublons, on renvoie une liste sans
+ doublons (et qui reviendra triée)
+ """
+
+ """ Au cas ou ca ne serait pas deja un vecteur numpy """ # Exemple
+ liste=N.sort(N.array(liste,copy=0)) # [1, 2, 2, 3, 3, 4, 5]
+ liste_diff=liste[1:]-liste[:-1] # [1, 0, 1, 0, 1, 1]
+ liste_temp=N.nonzero(liste_diff) # [0, 2, 4, 5]
+ liste_indice=N.zeros(liste_temp.shape[0]+1)
+ liste_indice[0]=0
+ liste_indice[1:]=liste_temp+1 # [0, 1, 3, 5, 6]
+ liste2=N.take(liste,liste_indice) # [1, 2, 3, 4, 5]
+ return liste2
+
+
+# ============================================================================ #
+
+
+class CONNEC(UserList.UserList) :
+ """
+ Connectivite : sequence mutable de sequences de noeuds
+ Pour l'instant, on s'appuie sur une liste de liste.
+ """
+
+ def __init__(self,nma) :
+
+ UserList.UserList.__init__(self,[None]*nma)
+
+
+ def Index(self) :
+ """
+ Retourne la connectivite sous forme de deux vecteurs numpy :
+ ind -> tableau des index (y compris le n+1 eme)
+ noe -> tableau des listes de noeuds
+ """
+
+ # Dimension des deux vecteurs
+ nma = len(self)
+ ltot = N.reduce(lambda x,y : x+len(y), self,0)
+ ind = N.zeros(nma+1,Int)
+ noe = N.zeros(ltot,Int)
+
+ # Construction des vecteurs
+ ind[0] = 0
+ for i in range(nma) :
+ m = self[i]
+ ind[i+1] = ind[i] + len(m)
+ noe[ind[i]:ind[i+1]] = N.array(m)
+
+ return ind,noe
+
+
+# ============================================================================ #
+# ============================================================================ #
+
+class MAIL_PY :
+
+ """
+ SD PYTHON MAILLAGE
+ La numeration est 0..N-1 pour les noeuds et 0..M-1 pour les mailles
+ """
+
+ def __init__(self,nno=0,nma=0) :
+
+ self.cn = N.zeros((nno,3),N.Float)
+ self.tm = N.zeros(nma,N.Int)
+ self.co = CONNEC(nma)
+ self.gma = {}
+ self.gno = {}
+ self.indice_noeuds = N.arange(nno)
+
+ self.correspondance_noeuds = []
+ self.correspondance_mailles = []
+
+ try:
+ import aster
+ nom_mailles = (None,) + string.strip(aster.getvectjev('&CATA.TM.NOMTM'))
+ except:
+ nom_mailles = (None,
+ 'POI1', 'SEG2', 'SEG22', 'SEG3', 'SEG33', 'SEG4', 'TRIA3',
+ 'TRIA33', 'TRIA6', 'TRIA66', 'TRIA7', 'QUAD4', 'QUAD44', 'QUAD8',
+ 'QUAD88', 'QUAD9', 'QUAD99', 'TETRA4', 'TETRA10','PENTA6', 'PENTA15',
+ 'PYRAM5', 'PYRAM13','HEXA8', 'HEXA20', 'HEXA27', 'TR3QU4', 'QU4TR3',
+ 'TR6TR3', 'TR3TR6', 'TR6QU4', 'QU4TR6', 'TR6QU8', 'QU8TR6', 'TR6QU9',
+ 'QU9TR6', 'QU8TR3', 'TR3QU8', 'QU8QU4', 'QU4QU8', 'QU8QU9', 'QU9QU8',
+ 'QU9QU4', 'QU4QU9', 'QU9TR3', 'TR3QU9', 'SEG32', 'SEG23' )
+
+ dic_mailles = {}
+ for i in range(len(nom_mailles)) :
+ dic_mailles[nom_mailles[i]] = i
+
+ self.nom = nom_mailles
+ self.dic = dic_mailles
+
+ try:
+ psyco.bind(self.FromAster)
+ except: pass
+
+
+# -------------------------------------------------------------
+
+ def get_connexite(self, nom, nma):
+ co=CONNEC(nma)
+ dico_connexite = aster.getcolljev(nom)
+ for element in dico_connexite.keys() :
+ co[int(element)-1] = (N.array(dico_connexite[element])-1)
+ return co
+
+
+ def get_coordonnees_noeuds(self, nom, nombre_noeuds):
+ coordonnees_noeuds = aster.getvectjev(nom)
+ coordonnees_noeuds = N.array(coordonnees_noeuds)
+ coordonnees_noeuds.shape = (nombre_noeuds,3)
+ cn = coordonnees_noeuds
+ return cn
+
+
+# -------------------------------------------------------------
+
+ def FromAster(self,nom) :
+
+ # On accepte le concept Aster ou bien la chaine texte de son nom
+ if type(nom)!=types.StringType:
+ nom_maillage = nom.nom
+ else:
+ nom_maillage = nom
+
+ nom_maillage=string.ljust(nom_maillage,8)
+
+ # recuperation de la taille
+ self.dime_maillage = aster.getvectjev(nom_maillage+'.DIME')
+ nombre_noeuds = self.dime_maillage[0]
+ nombre_mailles = self.dime_maillage[2]
+
+ # coordonnees des noeuds
+ self.cn = self.get_coordonnees_noeuds(nom_maillage+'.COORDO .VALE', nombre_noeuds)
+
+ # type de maille
+ self.tm = N.array(aster.getvectjev(nom_maillage+'.TYPMAIL'))
+
+ # connexite
+ self.co = self.get_connexite(nom_maillage+'.CONNEX', nombre_mailles)
+
+ self.indice_noeuds=N.arange(nombre_noeuds)
+
+ # groupe de noeuds
+ Lno_groupno_tot = aster.getcolljev(nom_maillage+'.GROUPENO')
+
+ Lno_groupno={}
+ for key in Lno_groupno_tot :
+ # Tolerance car parfois defi_group crée des groupes nuls à clé entiere
+ try:
+ Lno_groupno[key.strip()]=N.array(Lno_groupno_tot[key])-1
+ except: pass
+ self.gno=Lno_groupno
+
+ # groupe de mailles
+ Lma_groupma_tot = aster.getcolljev(nom_maillage+'.GROUPEMA')
+ Lma_groupma={}
+ for key in Lma_groupma_tot :
+ Lma_groupma[key.strip()]=N.array(Lma_groupma_tot[key])-1
+ self.gma=Lma_groupma
+
+ del(Lma_groupma_tot)
+
+ # listes de correspondance entre Aster et Mail-Py
+ self.correspondance_noeuds = aster.getvectjev(nom_maillage+'.NOMNOE')
+ self.correspondance_mailles = aster.getvectjev(nom_maillage+'.NOMMAI')
+
+
+# -------------------------------------------------------------
+
+ def ToAster(self) :
+
+ try:
+ INFO_EXEC_ASTER = self.jdc.get_cmd('INFO_EXEC_ASTER')
+ DETRUIRE = self.jdc.get_cmd('DETRUIRE')
+ LIRE_MAILLAGE = self.jdc.get_cmd('LIRE_MAILLAGE')
+ except:
+ try:
+ from Cata.cata import INFO_EXEC_ASTER
+ from Cata.cata import DETRUIRE
+ from Cata.cata import LIRE_MAILLAGE
+ except:
+ print "\n\nERREUR : il faut lancer ce programme depuis Aster pour pouvoir \ngénérer un maillage Aster.\n\n"
+ sys.exit(1)
+
+ # Recuperation d'une unité logique libre
+ _UL=INFO_EXEC_ASTER(LISTE_INFO='UNITE_LIBRE')
+ ul1=_UL['UNITE_LIBRE',1]
+ DETRUIRE(CONCEPT=(_F(NOM=_UL),), INFO=1)
+
+ # Sauvegarde du maillage dans un fichier .mail
+ fichier = 'fort.'+str(ul1)
+ f = open(fichier, 'w')
+ f.write( self.Voir_Mail() )
+ f.close()
+
+ # Récupération des concepts Aster présents dans la base
+# Note AA : ne marche pas encore bien :
+# Limitation : ne gere qu'un seul maillage genere par ToAster
+#
+# nom='_MSH'
+# try:
+# self.jdc_recupere=CONTEXT.get_current_step()
+# t_maillage=[]
+# for i in self.jdc_recupere.sds_dict.keys( ):
+# if self.jdc_recupere.sds_dict[i].__class__.__name__ == 'maillage_sdaster':
+# if (mail[0:4]==nom):
+# t_maillage.append( i )
+#
+# num=len(_lst)+1
+# except:
+# num=0
+
+ # Lecture de ce fichier .mail pour l'injecter dans l'espace Aster
+ _SMESH_ = LIRE_MAILLAGE(UNITE = ul1)
+
+ return(_SMESH_)
+
+# -------------------------------------------------------------
+
+ def __str__(self) :
+ return self.Voir_Mail()
+
+# -------------------------------------------------------------
+
+ def Voir_Mail(self) :
+ """
+ Impression au format ASTER
+ """
+
+ l = []
+ l.append('TITRE')
+ l.append('% CLASSE PY_MAIL -> MAIL')
+ l.append('FINSF')
+ l.append('%')
+
+ (nno,ndim) = self.cn.shape
+
+
+ # Coordonnees des noeuds
+ l.append('COOR_3D')
+
+ # Si le maillage initial ne provient pas d'Aster
+ if len(self.correspondance_noeuds) == 0:
+ for i in range(nno) :
+ ch = 'N'+repr(i)+' '+`self.cn[i,0]` + ' ' + `self.cn[i,1]` + ' ' + `self.cn[i,2]`
+# ch = 'N'+repr(i+1)+' '+`self.cn[i,0]` + ' ' + `self.cn[i,1]` + ' ' + `self.cn[i,2]`
+ l.append(ch)
+
+ # Si le maillage initial provient d'Aster
+ else:
+ for i in range(nno) :
+ ch = self.correspondance_noeuds[i]+' '+`self.cn[i,0]` + ' ' + `self.cn[i,1]` + ' ' + `self.cn[i,2]`
+# ch = 'N'+repr(i+1)+' '+`self.cn[i,0]` + ' ' + `self.cn[i,1]` + ' ' + `self.cn[i,2]`
+ l.append(ch)
+
+ # Connectivité des mailles
+ ind = N.argsort(self.tm)
+ ty = 0
+
+ # Si le maillage initial ne provient pas d'Aster
+ if len(self.correspondance_mailles) == 0:
+ for m in ind :
+ if self.tm[m] <> ty :
+ l.append('FINSF') ; l.append('%')
+ ty = self.tm[m]
+ l.append(self.nom[ty])
+ ch = 'M'+`m`+' '
+# ch = 'M'+`m+1`+' '
+ for n in self.co[m] :
+ ch = ch + 'N'+`n` + ' '
+# ch = ch + 'N'+`n+1` + ' '
+ l.append(ch)
+
+ # Si le maillage initial provient d'Aster
+ else:
+ for m in ind :
+ if self.tm[m] <> ty :
+ l.append('FINSF') ; l.append('%')
+ ty = self.tm[m]
+ l.append(self.nom[ty])
+ ch = self.correspondance_mailles[m]+' '
+# ch = 'M'+`m+1`+' '
+ for n in self.co[m] :
+ ch = ch + self.correspondance_noeuds[n] + ' '
+# ch = ch + 'N'+`n+1` + ' '
+ l.append(ch)
+
+ l.append('FINSF') ; l.append('%')
+
+
+ # Group_ma et Group_no
+ entete = ['GROUP_MA','GROUP_NO']
+ d_gp = [self.gma,self.gno]
+ pref = ['M','N']
+
+ # Si le maillage initial ne provient pas d'Aster
+ if len(self.correspondance_mailles) == 0:
+
+ for (d_gp,entete,prefixe) in [(self.gma,'GROUP_MA','M'),(self.gno,'GROUP_NO','N')] :
+ for gp in d_gp :
+ if len(d_gp[gp])>0: # On ne prend en compte que les group_* non vides
+ l.append(entete)
+ l.append(' ' + gp)
+ ch = ' '
+ i=0
+ for o in d_gp[gp]:
+ i+=1 # on ne met que 8 mailles sur une meme ligne
+ if (len(ch) > 60 or i>7):
+ l.append(ch)
+ ch = ' '
+ i=0
+ ch = ch + prefixe + `o` + ' '
+ l.append(ch)
+ l.append('FINSF') ; l.append('%')
+
+ # Si le maillage initial provient d'Aster
+ else:
+
+ for (d_gp,entete,prefixe) in [(self.gma,'GROUP_MA','M'),(self.gno,'GROUP_NO','N')] :
+ for gp in d_gp :
+ if len(d_gp[gp])>0: # On ne prend en compte que les group_* non vides
+ l.append(entete)
+ l.append(' ' + gp)
+ ch = ' '
+ i=0
+ for o in d_gp[gp]:
+ i+=1 # on ne met que 8 mailles sur une meme ligne
+ if (len(ch) > 60 or i>7):
+ l.append(ch)
+ ch = ' '
+ i=0
+# ch = ch + prefixe + `o` + ' '
+ if prefixe=='M':
+ ch = ch + self.correspondance_mailles[o] + ' '
+ else:
+ ch = ch + self.correspondance_noeuds[o] + ' '
+ l.append(ch)
+ l.append('FINSF') ; l.append('%')
+
+ # Fin
+ l.append('FIN\n')
+ return string.join(l,'\n')
+
+
+
+
+# ============================================================================ #
+# ============================================================================ #
+
+class PARTITION:
+
+ def __init__(self, jdc=None ,nb=0):
+
+ self.jdc = jdc
+
+ self.fichier_out = ''
+ self.liste_mailles = N.array( [] )
+ self.liste_sd = N.array( [] )
+ self.liste_mailles_bord = []
+ self.liste_sd_bord = []
+
+ self.MAILLAGE_Python = None
+
+ self.RELATIONS = { 'C_plus' : None,
+ 'C_moins': None,
+ 'nr': 0 }
+
+ self.ASTER = { 'MAILLAGE' : None,
+ 'MODELE' : None,
+ 'GROUP_MA' : None,
+ 'GROUP_MA_BORD' : None,
+ 'DICO_SD_MAILLES' : None,
+ }
+
+ self.OPTIONS = { 'NB_PART' : '',
+ 'ALGO' : '',
+ 'INFO' : '',
+ 'rep_metis' : aster.repout(),
+ 'exe_metis' : aster.repout() + 'pmetis',
+ 'fichier_in' : 'fort.66',
+ 'fichier_out' : 'fort.68',
+ 'elimine_bords': 'OUI',
+ }
+
+ self.Creation_Dico_Correspondance_Type_Maille()
+
+
+
+# ---------------------------------------------------------------------------- #
+
+ def __str__(self) :
+ """
+ Impression du contenu de la partition
+ """
+ l = []
+ l.append( 'Contenu de la partition :' )
+ l.append( '-------------------------' )
+ try: l.append( '- Maillage : ' + str(self.ASTER['MAILLAGE'].nom) )
+ except: pass
+ try: l.append( '- Modele : ' + str(self.ASTER['MODELE'].nom) )
+ except: pass
+ l.append( '- Nb part : ' + str(self.OPTIONS['NB_PART']) )
+ l.append( '- Niveau INFO : ' + str(self.OPTIONS['INFO']) )
+ l.append( '- Liste group_ma : ' + str(self.ASTER['GROUP_MA']) )
+
+ return string.join(l,'\n')
+
+# ---------------------------------------------------------------------------- #
+
+ def Partitionne_Aster(self, MAILLAGE, NB_PART, MODELE=None, METHODE=None, LOGICIEL=None, INFO=1):
+
+ self.t00 = time.clock()
+
+ self.OPTIONS['INFO'] = INFO
+
+ if MODELE:
+ # Recuperation de la liste des mailles à perndre en compte
+ self.ASTER['MODELE'] = MODELE
+ self.ASTER['MAILLAGE'] = MAILLAGE
+ _LST_MA = self.Modele_to_Liste_Mailles(MODELE)
+
+ elif MAILLAGE:
+ self.ASTER['MAILLAGE'] = MAILLAGE
+ _LST_MA = None
+
+
+ # Creation du maillage Python correspondant au maillage Aster
+ MAILLAGE_Python = MAIL_PY()
+ MAILLAGE_Python.FromAster(MAILLAGE.nom)
+
+ # Partitionne le maillage Python avec la liste de mailles _LST_MA
+ self.Partitionne_Maillage(MAILLAGE_Python, NB_PART, MAILLE=_LST_MA, METHODE=METHODE, LOGICIEL=LOGICIEL, INFO=INFO)
+
+ return
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Partitionne_Maillage(self, MAILLAGE_Python, NB_PART, MAILLE=None, METHODE=None, LOGICIEL=None, INFO=1):
+
+ self.t00 = time.clock()
+
+ if METHODE:
+ self.OPTIONS['exe_metis'] = aster.repout() + string.lower(METHODE)
+ elif LOGICIEL:
+ self.OPTIONS['exe_metis'] = LOGICIEL
+
+ self.OPTIONS['NB_PART'] = NB_PART
+ self.OPTIONS['INFO'] = INFO
+ self.MAILLAGE_Python = MAILLAGE_Python
+
+ exe_metis = self.OPTIONS['exe_metis']
+ f_metis = self.OPTIONS['fichier_in']
+ fw_metis = self.OPTIONS['fichier_out']
+
+ _LST_MA = MAILLE
+
+ # On initialise la connectivité et la connectivité inverse des aretes
+ self.MAILLAGE_Python.ca = {}
+ self.MAILLAGE_Python.cia = {}
+
+ _DIM = self.MAILLAGE_Python.dime_maillage[5]
+ _LST_TMA = self.MAILLAGE_Python.tm
+
+ if self.OPTIONS['INFO']>=5:
+ print 'cn=', self.MAILLAGE_Python.cn
+ print 'tm=', self.MAILLAGE_Python.tm
+ print 'co=', self.MAILLAGE_Python.co
+ print 'gma=', self.MAILLAGE_Python.gma
+ print 'gno=', self.MAILLAGE_Python.gno
+ print 'dim=', self.MAILLAGE_Python.dime_maillage
+ if self.OPTIONS['INFO']>=5: print '_LST_MA=', _LST_MA
+
+
+ # Elimination des mailles de bords
+ if self.OPTIONS['elimine_bords']!='NON':
+
+ # Liste des mailles à prendre en compte : dimension _DIM
+ _D_DIM_MAILLES = self.Creation_Listes_Mailles_Par_Dim(self.MAILLAGE_Python.tm, _LST_MA=_LST_MA)
+
+ # Connectivité et connectivité inverse sur les bords
+ self.Connectivite_Aretes()
+
+ self.liste_mailles = _D_DIM_MAILLES[ _DIM ]
+
+ # Pour prendre en compte des mélanges d'elements de dimension differente
+ _LST, _LST_BD = self.Elimination_Mailles_de_bords(MAILLAGE_Python, _D_DIM_MAILLES, _DIM)
+ self.liste_mailles = N.concatenate( (self.liste_mailles,N.array(_LST)) )
+
+ if self.OPTIONS['INFO']>=5:
+ print '_LST_BD=',_LST_BD
+ print '_LST=',_LST
+
+ else:
+ self.liste_mailles = _LST_MA
+
+
+ # Restriction des connectivités aux mailles à prendre en compte
+ self.Connectivite_Aretes(OPTION='all', _LST_OK=self.liste_mailles)
+
+ # Creation de l'arbre de connectivité des bords
+ self.Creation_Graphe()
+
+ # Reduction de l'arbre de connectivité des bords
+ _nb = self.Reduction_Graphe(_DIM)
+
+ # Ecriture du fichier pour Metis/Chaco/Jostle
+ _D_CORRES = self.Ecrire_Graphe(f_metis, _nb)
+
+ # Lancement de metis sur le fichier fort.UL (production de fort.UL.part.N)
+ txt = exe_metis + ' ' + f_metis + ' ' + str(NB_PART)
+ print 'Commande : ',txt
+ os.system( txt )
+
+ # Lecture du fichier resultant de Metis
+ self.fichier_out = f_metis + '.part.' + str(NB_PART)
+ self.liste_sd = self.Lecture_fichier_sdd(self.fichier_out, self.liste_mailles)
+
+ # Traitement des mailles de bords (on les reinjecte dans un SD)
+ if self.OPTIONS['elimine_bords']!='NON':
+ self.Affectation_Mailles_de_bords(_LST_BD, _DIM)
+
+ t1 = time.clock()
+ print "--- FIN PARTITIONNEMENT : ", t1 - self.t00
+
+ return
+
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Creation_Dico_Correspondance_Type_Maille(self):
+
+ # TYPE_ELEM : CF. &CATA.TM
+ # 1 - >POI1 <>SEG2 <>SEG22 <>SEG3 <>SEG33 <>SEG4 <>TRIA3 <
+ # 8 - >TRIA33 <>TRIA6 <>TRIA66 <>TRIA7 <>QUAD4 <>QUAD44 <>QUAD8 <
+ # 15 - >QUAD88 <>QUAD9 <>QUAD99 <>TETRA4 <>TETRA10 <>PENTA6 <>PENTA15 <
+ # 22 - >PYRAM5 <>PYRAM13 <>HEXA8 <>HEXA20 <>HEXA27 <>TR3QU4 <>QU4TR3 <
+ # 29 - >TR6TR3 <>TR3TR6 <>TR6QU4 <>QU4TR6 <>TR6QU8 <>QU8TR6 <>TR6QU9 <
+ # 36 - >QU9TR6 <>QU8TR3 <>TR3QU8 <>QU8QU4 <>QU4QU8 <>QU8QU9 <>QU9QU8 <
+ # 43 - >QU9QU4 <>QU4QU9 <>QU9TR3 <>TR3QU9 <>SEG32 <>SEG23 <
+
+ # Creation du dictionnaire des correspondance type_maille -> liste des aretes
+ maille2aretes={}
+ # POI
+ maille2aretes[1] = [ ]
+ # SEG
+ maille2aretes[2] = [ (0,1) ]
+ maille2aretes[3] = maille2aretes[4] = maille2aretes[5] = maille2aretes[6] = maille2aretes[2]
+ # TRIA
+ maille2aretes[7] = [ (0,1),(1,2),(0,2) ]
+ maille2aretes[8] = maille2aretes[9] = maille2aretes[10] = maille2aretes[11] = maille2aretes[7]
+ # QUAD
+ maille2aretes[12] = [ (0,1),(1,2),(2,3),(0,3) ]
+ maille2aretes[13] = maille2aretes[14] = maille2aretes[15] = maille2aretes[16] = maille2aretes[17] = maille2aretes[12]
+ # TETRA
+ maille2aretes[18] = [ (0,1,2),(0,1,3),(0,2,3),(1,3,2) ]
+ maille2aretes[19] = maille2aretes[18]
+ # PENTA
+ maille2aretes[20] = [ (0,1,2),(3,4,5),(0,2,5,3),(0,1,4,3),(2,1,4,5) ]
+ maille2aretes[21] = maille2aretes[20]
+ # PYRAM
+ maille2aretes[22] = [ (0,1,4),(1,2,4),(2,3,4),(3,0,4),(0,1,2,3) ]
+ maille2aretes[23] = maille2aretes[22]
+ # HEXA
+ maille2aretes[24] = [ (0,1,2,3), (4,5,6,7), (1,2,6,5), (2,3,7,6), (7,4,0,3), (4,5,1,0) ]
+ maille2aretes[25] = maille2aretes[26] = maille2aretes[24]
+
+
+ # dictionnaire de correspondance entre type_maille -> nb noeud (maille linéaire)
+ maille2nb={}
+ # POI
+ maille2nb[1] = 1
+ # SEG
+ maille2nb[2] = 2
+ maille2nb[3] = maille2nb[4] = maille2nb[5] = maille2nb[6] = maille2nb[2]
+ # TRIA
+ maille2nb[7] = 3
+ maille2nb[8] = maille2nb[9] = maille2nb[10] = maille2nb[11] = maille2nb[7]
+ # QUAD
+ maille2nb[12] = 4
+ maille2nb[13] = maille2nb[14] = maille2nb[15] = maille2nb[16] = maille2nb[17] = maille2nb[12]
+ # TETRA
+ maille2nb[18] = 4
+ maille2nb[19] = maille2nb[18]
+ # PENTA
+ maille2nb[20] = 5
+ maille2nb[21] = maille2nb[20]
+ # PYRAM
+ maille2nb[22] = 5
+ maille2nb[23] = maille2nb[22]
+ # HEXA
+ maille2nb[24] = 6
+ maille2nb[25] = maille2nb[26] = maille2nb[24]
+
+
+ # dictionnaire de correspondance entre type_maille -> dimension
+ maille2dim = {}
+ # POI
+ maille2dim[1] = 0
+ # SEG
+ maille2dim[2] = 1
+ maille2dim[3] = maille2dim[4] = maille2dim[5] = maille2dim[6] = maille2dim[2]
+ # TRIA
+ maille2dim[7] = 2
+ maille2dim[8] = maille2dim[9] = maille2dim[10] = maille2dim[11] = maille2dim[7]
+ # QUAD
+ maille2dim[12] = 2
+ maille2dim[13] = maille2dim[14] = maille2dim[15] = maille2dim[16] = maille2dim[17] = maille2dim[12]
+ # TETRA
+ maille2dim[18] = 3
+ maille2dim[19] = maille2dim[18]
+ # PENTA
+ maille2dim[20] = 3
+ maille2dim[21] = maille2dim[20]
+ # PYRAM
+ maille2dim[22] = 3
+ maille2dim[23] = maille2dim[22]
+ # HEXA
+ maille2dim[24] = 3
+ maille2dim[25] = maille2dim[26] = maille2dim[24]
+
+ # On stocke les dictionnaires
+ self.maille2aretes = maille2aretes
+ self.maille2nb = maille2nb
+ self.maille2dim = maille2dim
+
+ return
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Modele_to_Liste_Mailles(self, MODELE):
+
+ nommod = string.ljust(MODELE.nom,8)
+ _DIC_MA = aster.getcolljev(nommod.ljust(8)+'.MODELE .LIEL')
+
+ # Creation de la liste des mailles
+ ll = []
+ for type_maille in _DIC_MA.keys():
+ ll.extend( _DIC_MA[type_maille][0:-1] )
+ _LST_MA = N.array( ll ) - 1
+
+ if self.OPTIONS['INFO']>=5:
+ print '\n# ----- MODELE ----- #\n'
+ print '_LST_MA=',len(_LST_MA),_LST_MA
+ print '_DIC_MA=',len(_DIC_MA),_DIC_MA
+
+ return _LST_MA
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Creation_Listes_Mailles_Par_Dim(self, _LST_TMA, _LST_MA=None):
+
+ t0 = time.clock()
+
+ # Si _LST_MA est renseigné on extrait la liste des TMA correspondante aux mailles de _LST_MA
+ if _LST_MA != None:
+ _LST_TMA = N.take(_LST_TMA,_LST_MA)
+ else:
+ _LST_MA = N.arange(len(_LST_TMA))
+
+ _D_DIM_MAILLES = {}
+
+ # Liste des mailles 3D (type maille de 18 à 26)
+ _lst = N.where( _LST_TMA>=18, -3, _LST_TMA )
+ _tmp = N.where( _lst==-3, -1, 0 )
+# _D_DIM_MAILLES[3] = N.nonzero( _tmp )
+ _D_DIM_MAILLES[3] = N.take(_LST_MA, N.nonzero( _tmp ) )
+
+ # Liste des mailles 2D (type maille de 7 à 17)
+ _lst = N.where( _lst>=7, -2, _lst )
+ _tmp = N.where( _lst==-2, -1, 0 )
+ _D_DIM_MAILLES[2] = N.take(_LST_MA, N.nonzero( _tmp ) )
+
+ # Liste des mailles 1D (type maille de 2 à 6)
+ _lst = N.where( _lst>=2, -1, _lst )
+ _tmp = N.where( _lst==-1, -1, 0 )
+ _D_DIM_MAILLES[1] = N.take(_LST_MA, N.nonzero( _tmp ) )
+
+ # Liste des mailles 0D (type maille 1)
+ _lst = N.where( _lst>=1, -4, _lst )
+ _tmp = N.where( _lst==-4, -1, 0 )
+ _D_DIM_MAILLES[0] = N.take(_LST_MA, N.nonzero( _tmp ) )
+
+
+ if self.OPTIONS['INFO']>=5:
+ for i in _D_DIM_MAILLES.keys():
+ print "-----------------"
+ print 'Dim:',i, _D_DIM_MAILLES[i]
+ print "-----------------"
+
+ print "--- FIN Creation des listes de mailles par Dim : ", time.clock() - t0
+
+ return _D_DIM_MAILLES
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Connectivite_Aretes(self, OPTION=None, _LST_OK=None):
+
+ t0 = time.clock()
+
+ # Si _LST_OK n'est pas renseigné on prend toutes les mailles
+ if not _LST_OK: _LST_OK = N.arange(len(self.MAILLAGE_Python.tm))
+
+ if self.OPTIONS['INFO']>=5: print '_LST_OK (ca)=',_LST_OK
+
+ maille2aretes = self.maille2aretes
+
+ # Creation de la :
+ # - connectivite des aretes (self.MAILLAGE_Python.ca) : m1 -> [ (a1, a2), .. ]
+ # - connectivite inverse des aretes (self.MAILLAGE_Python.cia) : (a1, a2) -> [ m1, m2, ... ]
+
+ self.MAILLAGE_Python.ca = {}
+ self.MAILLAGE_Python.cia = {}
+
+ for n in _LST_OK:
+
+ n1 = self.MAILLAGE_Python.tm[n]
+
+ l_aretes = maille2aretes[n1] # liste des aretes de la maille n
+ l_noeuds = self.MAILLAGE_Python.co[n] # liste des noeuds de la maille n
+
+ for arete in l_aretes:
+ ll = []
+ for i in arete:
+ ll.append( l_noeuds[i] )
+ ll.sort()
+ ll = tuple(ll)
+
+ # Table de connectivité des aretes
+ if OPTION:
+ if not self.MAILLAGE_Python.ca.has_key(n): self.MAILLAGE_Python.ca[n]=[]
+ self.MAILLAGE_Python.ca[n].append(ll)
+# try:
+# self.MAILLAGE_Python.ca[n].append(ll)
+# except KeyError:
+# self.MAILLAGE_Python.ca[n]=[ll]
+
+ # Table de connectivité inverse des aretes
+ if not self.MAILLAGE_Python.cia.has_key(ll): self.MAILLAGE_Python.cia[ll]=[]
+ self.MAILLAGE_Python.cia[ll].append(n)
+# try:
+# self.MAILLAGE_Python.cia[ll].append(n)
+# except KeyError:
+# self.MAILLAGE_Python.cia[ll]=[n]
+
+
+ if self.OPTIONS['INFO']>=5:
+ for k in self.MAILLAGE_Python.cia.keys():
+ print 'cia:',k, ' ', self.MAILLAGE_Python.cia[k]
+ if OPTION:
+ for k in self.MAILLAGE_Python.ca.keys():
+ print 'ca: ',k, ' ', self.MAILLAGE_Python.ca[k]
+
+
+ print "--- FIN Creation de la connectivite simple et inverse des aretes : ", time.clock() - t0
+
+ return
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Elimination_Mailles_de_bords(self, MAILLAGE_Python, _D_DIM_MAILLES, _DIM):
+ """
+ Extraction des mailles de bords (mailles incluses dans un bord d une autre maille)
+ """
+
+ t0 = time.clock()
+
+ _LST_TMA = self.MAILLAGE_Python.tm
+
+ if self.OPTIONS['INFO']>=5:
+ MAILLAGE = self.ASTER['MAILLAGE']
+ nommail = string.ljust(MAILLAGE.nom,8)
+ _LST_MAI = aster.getvectjev(nommail.ljust(8)+'.NOMMAI')
+
+ # Le dico maille2nb donne le nombre de noeuds definissant un bord (lineaire)
+ maille2nb = self.maille2nb
+
+
+ # construction des listes des mailles de dim N-1 :
+ # _LST_OK : Mailles de dim N-i qui ne sont pas un bord des mailles de dim N
+ # _LST_BD : Mailles de dim N-i qui sont un bord
+ #
+ if self.OPTIONS['INFO']>=5: print '\n\nElimination des mailles de bord de DIM', _DIM - 1
+
+ _LST4 = _D_DIM_MAILLES[ _DIM - 1 ]
+ _LST_IND = N.arange( len(_LST4) ) + 1 # on ajoute 1 pour eviter le premier 0 dans les test nonzero plus bas
+
+ if self.OPTIONS['INFO']>=5: print ' Mailles concernées=',_LST4
+
+ i=0
+ for m in _LST4:
+ if self.OPTIONS['INFO']>=5: print '\n Maille de dim N-1:',m, ' Aster:',string.strip(_LST_MAI[m]), ' TMA:',self.MAILLAGE_Python.tm[m], ' CO:',self.MAILLAGE_Python.co[m], '(noeuds de cette maille)'
+ nb = maille2nb[ self.MAILLAGE_Python.tm[m] ]
+ ll = self.MAILLAGE_Python.co[m][0:nb]
+ ll = N.sort(ll)
+ ll = ll.tolist()
+ ll = tuple(ll)
+ if self.OPTIONS['INFO']>=5: print ' Bord (lineaire)', ll, nb
+
+ try:
+ if self.OPTIONS['INFO']>=5: print ' CIA=', self.MAILLAGE_Python.cia[ ll ], '(mailles de dim N qui ont cette maille pour bord)'
+ _tmp=[]
+ for maille in self.MAILLAGE_Python.cia[ ll ]:
+ if self.OPTIONS['INFO']>=5: print ' Maille N:', maille, 'Aster:', string.strip(_LST_MAI[maille]), ' TMA:', self.MAILLAGE_Python.tm[maille]
+# self.liste_mailles_bord.append(m)
+ except:
+ if self.OPTIONS['INFO']>=5: print ' Maille non-bord'
+ _LST_IND[i] = 0
+
+ i+=1
+
+ # Recuperation des mailles de bords et non-bords
+ _LST_BD = N.nonzero(_LST_IND)
+ _LST_BD = N.take(_LST4,_LST_BD)
+
+ _LST_OK = N.where( _LST_IND==0, 1 , 0 )
+ _LST_OK = N.nonzero(_LST_OK)
+ _LST_OK = N.take(_LST4,_LST_OK)
+
+ if self.OPTIONS['INFO']>=5: print '\nListe Maille de bords de DIM', _DIM - 1,' :',_LST_BD
+ if self.OPTIONS['INFO']>=5: print 'Liste Maille de DIM', _DIM - 1,'qui ne sont pas des bords :',_LST_OK
+
+ print "--- FIN Maille de bords de DIM",_DIM - 1, " : ", time.clock() - t0
+ t0 = time.clock()
+
+
+ # On cherche à marier les mailles de dimension N-2, N-3
+ # Peut etre lent car on utilise la connectivité ! Mais pour le moment on a rien d'autre.
+
+ _LST_BD0 = []
+ _LST_OK0 = []
+ _D_BD = {}
+ for d in range(_DIM-1):
+ _LST4 = _D_DIM_MAILLES[ d ]
+ if self.OPTIONS['INFO']>=5: print '\n\nElimination des mailles de bord de DIM', d
+ if self.OPTIONS['INFO']>=5: print ' Mailles concernées=',_LST4
+ for mai in _LST4:
+ if self.OPTIONS['INFO']>=5: print '\n Maille:', mai, ' Aster:',string.strip(_LST_MAI[mai]), ' TMA:',self.MAILLAGE_Python.tm[mai], ' CO:',self.MAILLAGE_Python.co[mai], '(noeuds de cette maille)'
+
+ nb = maille2nb[ self.MAILLAGE_Python.tm[mai] ]
+ ll = self.MAILLAGE_Python.co[mai][0:nb]
+ ll = N.sort(ll)
+ ll = ll.tolist()
+ _tmp = tuple(ll)
+# _tmp = self.MAILLAGE_Python.co[mai]
+
+ if self.OPTIONS['INFO']>=5: print ' Bord (lineaire):', _tmp, nb
+
+ ok=0
+ for arete in self.MAILLAGE_Python.cia:
+ _nb=0
+ for noe in _tmp:
+ if noe in arete: _nb+=1
+ if _nb == len(_tmp):
+ if self.OPTIONS['INFO']>=5: print ' Maille N+i:', self.MAILLAGE_Python.cia[arete], '- Arete:', arete
+ _LST_BD0.append( mai )
+ ok=1
+# if not _D_BD.has_key( mai ): _D_BD[mai] = []
+# _D_BD[mai].append( self.MAILLAGE_Python.cia[arete] )
+ break
+ if ok == 0:
+ _LST_OK0.append( mai )
+
+# print 'Mai:',mai, '_D_BD[mai]=',_D_BD[mai]
+
+
+ if self.OPTIONS['INFO']>=5: print '\nListe Maille de bords de DIM', d,' :',_LST_BD0
+ if self.OPTIONS['INFO']>=5: print 'Liste Maille de DIM', d,'qui ne sont pas des bords :',_LST_OK0
+
+
+ print '--- FIN Maille de bords de DIM', d, ' :',time.clock() - t0
+ t0 = time.clock()
+
+
+ _LST_OK = N.concatenate( (_LST_OK, N.array(_LST_OK0)) )
+ _LST_BD = N.concatenate( (_LST_BD, N.array(_LST_BD0)) )
+
+ if self.OPTIONS['INFO']>=5: print '\nTotal:\nListe Maille de bords=',_LST_BD
+ if self.OPTIONS['INFO']>=5: print 'Liste Maille non-bords=',_LST_OK,'\n'
+
+# print "--- FIN Maille de bords 3 : ", time.clock() - t0
+
+ return _LST_OK, _LST_BD
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Affectation_Mailles_de_bords(self, _LST_BD, _DIM):
+ """
+ Affectation a un SD des mailles de bords (mailles incluses dans un bord d une autre maille)
+ """
+
+ if self.OPTIONS['INFO']>=5:
+ print 'liste_mailles_bord=', self.liste_mailles_bord
+ print 'liste_sd_bord', self.liste_sd_bord
+ print '_LST_BD=',_LST_BD
+
+
+ MAILLAGE = self.ASTER['MAILLAGE']
+ _LST_TMA = self.MAILLAGE_Python.tm
+
+ if self.OPTIONS['INFO']>=5:
+ nommail = string.ljust(MAILLAGE.nom,8)
+ _LST_MAI = aster.getvectjev(nommail.ljust(8)+'.NOMMAI')
+
+ t0 = time.clock()
+
+ # Affectation des mailles de bords à chacun des SD
+
+ # Le dico maille2nb donne le nombre de noeuds definissant un bord (lineaire)
+ maille2nb = self.maille2nb
+
+ i = 0
+ for m in _LST_BD:
+ if self.OPTIONS['INFO']>=5: print '\n Maille de dim N-1:',m, ' Aster:',string.strip(_LST_MAI[m]), ' TMA:',self.MAILLAGE_Python.tm[m], ' CO:',self.MAILLAGE_Python.co[m], '(noeuds de cette maille)'
+ nb = maille2nb[ self.MAILLAGE_Python.tm[m] ]
+ ll = self.MAILLAGE_Python.co[m][0:nb]
+ ll = N.sort(ll)
+ ll = ll.tolist()
+ ll = tuple(ll)
+ if self.OPTIONS['INFO']>=5: print ' Bord (lineaire)', ll, nb
+
+ # Cas particulier des POI1 en 2D et 3D (ils ne peuvent etre des bords d'elements 2D ou 3D)
+ if ( (nb==1) and (_DIM>=2) ):
+ _tmp=[]
+ for arete in self.MAILLAGE_Python.cia.keys():
+ if ll[0] in arete:
+ for maille in self.MAILLAGE_Python.cia[ arete ]:
+ if self.OPTIONS['INFO']>=5: print ' Maille N+i:', maille, ' Aster:',string.strip(_LST_MAI[maille]), ' Arete:', arete
+ _tmp.append( self.liste_sd[maille] )
+
+ # Cas particulier des SEG en 3D (ils ne peuvent etre des bords d'elements 3D)
+ elif ( (nb==2) and (_DIM==3) ):
+ _tmp=[]
+ for arete in self.MAILLAGE_Python.cia.keys():
+ _nb=0
+ for noe in ll:
+ if noe in arete: _nb+=1
+ if _nb == len(ll):
+ for maille in self.MAILLAGE_Python.cia[arete]:
+ if self.OPTIONS['INFO']>=5: print ' Mailles N+i:', maille, ' Aster:',string.strip(_LST_MAI[maille]), ' Arete:', arete
+ _tmp.append( self.liste_sd[maille] )
+
+ # Autres mailles de bord
+ else:
+ if self.OPTIONS['INFO']>=5: print ' CIA=', self.MAILLAGE_Python.cia[ ll ], '(mailles de dim N qui ont cette maille pour bord)'
+ _tmp=[]
+ for maille in self.MAILLAGE_Python.cia[ ll ]:
+ if self.OPTIONS['INFO']>=5: print ' Maille N+i:', maille, 'Aster:', string.strip(_LST_MAI[maille]), ' SD:', self.liste_sd[maille], ' TMA:', self.MAILLAGE_Python.tm[maille]
+ _tmp.append( self.liste_sd[maille] )
+
+ # integre la maille au SD le plus faible (pour que des groupes de bords se retrouvent dans le meme SD)
+ _tmp.sort()
+ self.liste_mailles_bord.append(m)
+ self.liste_sd_bord.append( _tmp[0] )
+ i += 1
+ if self.OPTIONS['INFO']>=5: print ' ---> Maille:',m,'integree au SD:', _tmp[0]
+
+ if self.OPTIONS['INFO']>=5:
+ print '\n\nliste_mailles_bord=', self.liste_mailles_bord
+ print 'liste_sd_bord=', self.liste_sd_bord
+
+
+ print "--- FIN Affectation des mailles de bords : ", time.clock() - t0
+
+ return
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Creation_Graphe(self):
+
+ t0 = time.clock()
+
+ # Creation du graphe complet
+ self.GRAPH = {}
+
+ for mai in self.liste_mailles:
+ _ll=[]
+ for are in self.MAILLAGE_Python.ca[mai]:
+ _ll.extend( self.MAILLAGE_Python.cia[are] )
+ _mm = enleve_doublons_liste(_ll) # coute cher!
+ _tmp = _mm.tolist()
+ _tmp.remove(mai)
+ self.GRAPH[mai] = _tmp
+
+ if self.OPTIONS['INFO']>=5: print 'self.GRAPH['+str(mai)+']=', self.GRAPH[mai]
+
+ print "--- FIN Creation du graphe complet : ", time.clock() - t0
+
+ return
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Reduction_Graphe(self, _DIM):
+
+ t0 = time.clock()
+
+ # Elimination des connectivités à interface nulle
+ maille2dim = self.maille2dim
+ _lst2 = []
+ for mai in self.liste_mailles:
+ if self.OPTIONS['INFO']>=5: print '\nmai:', mai, 'co:', self.MAILLAGE_Python.co[mai], 'tm:', self.MAILLAGE_Python.tm[mai]
+ _DIM1 = maille2dim[ self.MAILLAGE_Python.tm[mai] ]
+ _tmp2 =[]
+ for mai2 in self.GRAPH[mai]:
+ if self.OPTIONS['INFO']>=5: print 'mai2:', mai2, 'co:', self.MAILLAGE_Python.co[mai2], 'tm:', self.MAILLAGE_Python.tm[mai2]
+ # calcule le nombre de noeuds communs aux deux mailles
+ _nb = 0
+ for noe in self.MAILLAGE_Python.co[mai2]:
+ if noe in self.MAILLAGE_Python.co[mai]: _nb += 1
+ _DIM2 = maille2dim[ self.MAILLAGE_Python.tm[mai2] ]
+ if _nb >= min(_DIM1, _DIM2): # le min permet de faire du collage 3D-coque par exemple
+ _tmp2.append( mai2 )
+ _tmp = [mai, mai2]
+ _tmp.sort()
+ _lst2.append(_tmp)
+ self.GRAPH[mai] = _tmp2
+
+ print "--- FIN Elimination des connectivités avec une interface nulle : ", time.clock() - t0
+ t0 = time.clock()
+
+
+ # Calcul du nombre d'aretes
+ # A voir : normalement il n'y a rien a faire car nb0 = 2*nb (a verifier...)
+ _lst2.sort()
+ _v = _lst2[0]
+ _nb = 1
+ for i in _lst2:
+ if i != _v:
+ _v = i
+ _nb += 1
+
+
+ if self.OPTIONS['INFO']>=5:
+ print '----------------------------------------------'
+ for mai in self.liste_mailles:
+ print 'self.GRAPH['+str(mai)+']=', self.GRAPH[mai]
+ print '----------------------------------------------'
+
+ return _nb
+
+
+# ------------------------------------------------------------------ #
+
+ def Ecrire_Graphe(self, f_metis, _nb):
+
+ t0 = time.clock()
+
+ # On doit renumeroter les mailles qui arrivent dans self.liste_mailles pour avoir 0... N-1
+ _D_CORRES = {}
+ for i in N.arange(len(self.liste_mailles)):
+ _D_CORRES[ self.liste_mailles[i] ] = i
+
+ # Ecriture du fichier fort.UL pour metis
+ fw = open(f_metis,'w')
+ fw.write( str(len(self.liste_mailles)) + ' ' + str(_nb) + '\n')
+ for l in self.liste_mailles:
+# try:
+ _tmp = []
+ for t in self.GRAPH[l]:
+ try:
+ t = _D_CORRES[t]
+ _tmp.append( str(t+1) ) # Necessaire car metis numerote de 1 à N
+ except:
+ print 'on oublie le bord:', t
+ fw.write( string.join(_tmp, ' ') + '\n' )
+# except:
+# print 'Probleme ecriture graphe! On continue..'
+ fw.close()
+
+ print "--- FIN Ecriture du fichier du graphe pour metis : ", time.clock() - t0
+
+ return _D_CORRES
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Lecture_fichier_sdd(self, fichier, _LST_OK):
+
+ t0 = time.clock()
+
+ # Lecture du fichier produit par metis (partie a optimiser)
+ try:
+ f = open( fichier, 'r' )
+ except:
+ print "\n\n ERREUR: le fichier est introuvable! Le partitionneur \n ne s'est probablement pas lancé.\n\n"
+ sys.exit(1)
+ else:
+ _tmp = []
+ for l in f.readlines():
+ _tmp.append( int(string.strip(l)) )
+ f.close()
+ _l_domaines = N.array(_tmp,copy=0)
+
+ # Pour garder le fichier metis
+ os.system( 'mv ' + fichier + ' REPE_OUT/' )
+
+ if self.OPTIONS['INFO']>=5: print '_l_domaines=',_l_domaines
+
+ print "--- FIN Lecture du fichier produit par metis : ", time.clock() - t0
+
+ return _l_domaines
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Creation_Group_ma_Python_par_SD(self, NOM='SD', NOM2='B'):
+
+ t0 = time.clock()
+
+ NB_PART = self.OPTIONS['NB_PART']
+
+ # Creation du dictionnaire des listes des mailles par SD
+ # d_gma : { num sd -> [ liste mailles ] }
+ d_gma = {}
+ for i in range(NB_PART):
+ d_gma[i] = []
+
+ i=0
+ for sdd in self.liste_sd:
+ d_gma[sdd].append( self.liste_mailles[i] )
+ i+=1
+
+
+ # Creation du dictionnaire des listes des mailles de bord par SD
+ # d_gma_bord : { num sd -> [ liste mailles ] }
+ d_gma_bord = {}
+ for i in range(NB_PART):
+ d_gma_bord[i] = []
+
+ i=0
+ for sdd in self.liste_sd_bord:
+ d_gma_bord[sdd].append( self.liste_mailles_bord[i] )
+ i+=1
+
+
+ # Generation des listes de noms de groupes
+ _l_sd = []
+ _l_bord = []
+ for i in range(NB_PART):
+ if d_gma[i] != []:
+ _l_sd.append( NOM + str(i) )
+ if d_gma_bord[i] != []:
+ _l_bord.append( NOM2 + str(i) )
+
+ # Stockage
+ self.ASTER['GROUP_MA'] = _l_sd
+ self.ASTER['GROUP_MA_BORD'] = _l_bord
+
+
+ # Creation des groupes de mailles dans le Maillage Python
+ for i in range(NB_PART):
+ self.MAILLAGE_Python.gma[NOM+str(i)] = d_gma[i]
+ self.MAILLAGE_Python.gma[NOM2+str(i)] = d_gma_bord[i]
+
+ print "--- FIN creation du dictionnaire des listes des mailles par SD ", time.clock() - t0
+
+ return
+
+
+# ---------------------------------------------------------------------------- #
+
+ def Creation_Group_ma_Aster_par_SD(self, NOM='SD', NOM2='B', INCLUSE='NON'):
+
+ t0 = time.clock()
+
+ MAILLAGE = self.ASTER['MAILLAGE']
+ NB_PART = self.OPTIONS['NB_PART']
+
+ nommail = string.ljust(MAILLAGE.nom,8)
+ _LST_MAI = aster.getvectjev(nommail.ljust(8)+'.NOMMAI')
+
+
+ # Creation du dictionnaire des listes des mailles par SD
+ # d_gma : { num sd -> [ liste mailles ] }
+ d_gma = {}
+ for i in range(NB_PART):
+ d_gma[i] = []
+
+ m=0
+ for sdd in self.liste_sd:
+ d_gma[sdd].append( string.strip(_LST_MAI[ self.liste_mailles[m] ]) ) # voir si le strip coute cher !
+ m += 1
+
+
+ # Creation du dictionnaire des listes des mailles de bord par SD
+ # d_gma_bord : { num sd -> [ liste mailles ] }
+ d_gma_bord = {}
+ for i in range(NB_PART):
+ d_gma_bord[i] = []
+
+ # On inclus directement les mailles de bords dans les SD
+ if INCLUSE=='OUI':
+ m=0
+ for sdd in self.liste_sd_bord:
+ d_gma[sdd].append( string.strip(_LST_MAI[ self.liste_mailles_bord[m] ]) ) # voir si le strip coute cher !
+ m+=1
+
+ else:
+ m=0
+ for sdd in self.liste_sd_bord:
+ d_gma_bord[sdd].append( string.strip(_LST_MAI[ self.liste_mailles_bord[m] ]) ) # voir si le strip coute cher !
+ m+=1
+
+
+ print "--- FIN creation du dictionnaire des listes des mailles par SD ", time.clock() - t0
+ t0 = time.clock()
+
+
+ # Creation et lancement de la commande DEFI_GROUP associée
+ try:
+ DEFI_GROUP = self.jdc.get_cmd('DEFI_GROUP')
+ except:
+ try:
+ from Cata.cata import DEFI_GROUP
+ except:
+ print "\n\nERREUR : il faut lancer ce programme depuis Aster pour pouvoir \ngénérer les groupes de mailles Aster.\n\n"
+ return
+
+ _tmp = []
+ _l_sd = []
+ _l_bord = []
+ for i in range(NB_PART):
+ if d_gma[i] != []:
+ _tmp.append( {'MAILLE': d_gma[i],'NOM': NOM + str(i)} )
+ _l_sd.append( NOM + str(i) )
+ if d_gma_bord[i] != []:
+ _tmp.append( {'MAILLE': d_gma_bord[i],'NOM': NOM2 + str(i)} )
+ _l_bord.append( NOM2 + str(i) )
+
+ motscle2= {'CREA_GROUP_MA': _tmp }
+
+ DEFI_GROUP( reuse=MAILLAGE,
+ MAILLAGE=MAILLAGE,
+ INFO=1,
+ **motscle2
+ ) ;
+
+ # Stockage
+ self.ASTER['DICO_SD_MAILLES'] = d_gma
+ self.ASTER['GROUP_MA'] = _l_sd
+ self.ASTER['GROUP_MA_BORD'] = _l_bord
+
+ print "--- FIN Creation et lancement de la commande DEFI_GROUP associée : ", time.clock() - t0
+
+ return
+
+# ---------------------------------------------------------------------------- #
--- /dev/null
+#@ MODIF sup_gmsh Utilitai DATE 10/05/2005 AUTEUR GJBHHEL E.LORENTZ
+# -*- coding: iso-8859-1 -*-
+# CONFIGURATION MANAGEMENT OF EDF VERSION
+# ======================================================================
+# COPYRIGHT (C) 1991 - 2004 EDF R&D WWW.CODE-ASTER.ORG
+# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
+# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
+# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
+# (AT YOUR OPTION) ANY LATER VERSION.
+#
+# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
+# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
+# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
+# GENERAL PUBLIC LICENSE FOR MORE DETAILS.
+#
+# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
+# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
+# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
+# ======================================================================
+
+import os.path, string, os, copy
+from Numeric import *
+
+try :
+ from Cata.cata import *
+ from Accas import _F
+except :
+ print 'Fonctionnalites Aster indisponibles'
+
+
+_CARAC = {
+ 'prec' : 1.E-8
+ }
+
+
+def Class_search(class_name, target_class) :
+
+ """
+ Check that class_name inherits from target_class
+ (run recursively through the inheritance lists)
+ """
+
+ if class_name == target_class : return 1
+
+ for cl in class_name.__bases__ :
+ if Class_search(cl, target_class) : return 1
+ return 0
+
+
+def Progress(L,**para) :
+
+ """
+ Compute the unknown parameters for a geometric progression :
+ r = ratio
+ N = number of elements
+ h = initial size
+
+ So that :
+ L = Sum(i=0:N-1, h(i)) where h(i+1) = h(i)*r, h(0)=h
+
+ Usage :
+ Progress(L,r=...,N=...) -> h
+ Progress(L,r=...,h=...) -> N
+ Progress(L,h=...,N=...) -> r
+
+ """
+
+ prec = 1.E-4
+
+ # Calcul de N
+ if 'N' not in para.keys() :
+ r = float(para['r'])
+ h = float(para['h'])
+ N = log(1+(r-1)*L/h)/log(r)
+ N = int(N+0.5)
+ return N
+
+ # Calcul de h
+ elif 'h' not in para.keys() :
+ r = float(para['r'])
+ N = int(para['N']+0.5)
+ h = L*(r-1)/(r**N-1)
+ return h
+
+ # Calcul de r
+ elif 'r' not in para.keys() :
+ h = float(para['h'])
+ N = int(para['N']+0.5)
+ a = L/h
+ if N > a :
+ x = 0
+ else :
+ x = a**(1./(N-1))
+
+ for i in xrange(100) :
+ res = x**N - a*x + a-1
+ if abs(res) < prec*(x-1)*a :
+ return x
+ dx = res/(a-N*x**(N-1))
+ x = x+dx
+
+ raise 'Solution failure'
+
+ else :
+ raise 'Unknown parameters'
+
+
+
+def Is_Geometric(object) :
+
+ """
+ return true if the object inherits of the Geometric class
+ """
+
+ return Class_search(object.__class__, Geometric)
+
+
+
+# --------------------------------------------------------------------------------------
+
+class Geometric :
+
+ """
+ GENERIC CLASS FOR GEOMETRICAL OBJECTS
+
+ private attribute
+ parameters : dictionnary of the attributes (except relation and parameters itself)
+ see __getattr__ and __setattr__
+
+
+ Attributes
+ num : index among gmsh objects
+ md : mesh descriptor
+ mesh : related mesh object
+ relation : model object in case of coincidence
+ type : type of the object (None, Point, Line, Circle, Surface, Volume)
+
+ Public methods
+ Is_point : return true is the object inherits of the Point class
+
+ Is_line : return true is the object inherits of the Line class
+
+ Is_surface : return true is the object inherits of the Surface class
+
+ Is_volume : return true is the object inherits of the Volume class
+
+ Base_class : return the name of the base class : Point, Line, Surface
+ or Volume
+
+ Is_same_dimension : return true is both objects are of the same dimension
+ (point, line, surface or volume)
+ in -> object to compare to self
+
+ Duplicate : duplicate an object and base its mesh_descriptor
+ on the mesh_descriptor of the model
+
+ Coincide : assert that an object is coincident with a model one
+ All the attributes are then automatically read from
+ the model object (see __setattr__ and __getattr__).
+ in -> model object
+
+ Private method
+
+ Root :
+ Provides the root object of an object, ie the object itself if there is no relation
+ or the deepest model in case of relation.
+
+ Geometric_coincide : check if a geometrical coincidence is possible
+ return information about the coincidence, false else.
+ in -> model object
+
+ Deep_coincide : proceed recursively to ensure coincidence of the relevant sub-objects
+ in -> model object
+ in -> correspond (information returned by Geometric_coincide)
+
+ __setattr__ : distinguish two sets of attributes
+ relation (to express a relation with a model object in case of coincidence)
+ all the other attributes which are stored in the dictionnary parameters
+ instead of the usual __dict__ if there is no relation (see Coincide)
+ and in the model object if there is a coincidence
+
+ __getattr__ : if the object is related (relation <> None) the attribute is read
+ in the model object. Else, it is read in the current object, actually
+ in the dictionnary parameters (see __setattr__)
+
+ Thanks to these two overloaded methods, the access to the attributes is usual if
+ there is no relation whereas the attributes of the model object are accessed
+ transparently if there is a relation.
+
+ __cmp__ :
+ The comparison of two objects involves possible coincidence. It is no more the object ids
+ that are compared but the object roots (.relation if any).
+
+ Gmsh : produce the source code for Gmsh
+ in -> mesh
+
+ Gmsh_send : send a line code to the gmsh interpreter
+ in -> line_code (string)
+
+ Intermediate_meshing : produce the source code for the intermediate objects
+ in -> mesh
+
+ Object meshing : produce the source code for the current object
+ var -> object number (modified if several objects are created)
+
+ """
+
+ def __init__(self) :
+ self.relation = None
+ self.parameters = {}
+ self.num = 0
+ self.md = Mesh_Descriptor()
+
+ types = {}
+ types[Geometric] = 'Geometric'
+ types[Point] = 'Point'
+ types[Line] = 'Line'
+ types[Circle] = 'Circle'
+ types[Surface] = 'Surface'
+ types[Volume] = 'Volume'
+ types[LineLoop] = 'LineLoop'
+ types[SurfaceLoop] = 'SurfaceLoop'
+ try :
+ self.type = types[self.__class__]
+ except KeyError :
+ raise 'Unknown object type'
+
+
+ def Is_point(self) :
+ return Class_search(self.__class__, Point)
+
+ def Is_line(self) :
+ return Class_search(self.__class__, Line)
+
+ def Is_surface(self) :
+ return Class_search(self.__class__, Surface)
+
+ def Is_volume(self) :
+ return Class_search(self.__class__, Volume)
+
+ def Base_class(self) :
+ if self.Is_volume() : return 'Volume'
+ if self.Is_surface() : return 'Surface'
+ if self.Is_line() : return 'Line'
+ if self.Is_point() : return 'Point'
+
+ def Is_same_dimension(self, obj) :
+
+ return (
+ (self.Is_point() and obj.Is_point() ) or
+ (self.Is_line() and obj.Is_line() ) or
+ (self.Is_surface() and obj.Is_surface() ) or
+ (self.Is_volume() and obj.Is_volume() )
+ )
+
+
+ def __setattr__(self, attr, value) :
+
+ if attr in ['relation','parameters'] :
+ self.__dict__[attr] = value
+ else :
+ if self.relation :
+ setattr(self.relation,attr,value)
+ else :
+ self.parameters[attr] = value
+
+
+ def __getattr__(self,attr) :
+
+ if self.relation :
+ return (getattr(self.relation,attr))
+ else :
+ if attr in self.parameters.keys() :
+ return self.parameters[attr]
+ else :
+ raise AttributeError,attr
+
+
+ def Root(self) :
+
+ o = self
+ while o.relation : o = o.relation
+ return o
+
+
+ def __cmp__(self,obj) :
+
+ if self.Root() is obj.Root() :
+ return 0
+ else :
+ return -1
+
+
+ def Geometric_coincide(self,obj) : return 0
+
+ def Deep_coincide(self,obj,correspond) : pass
+
+ def Coincide(self, obj) :
+
+ if self == obj : return # in that way recursive loops cannot exist
+
+ if self.relation : # the root is put in coincidence, not the object itself
+ self.Root().Coincide(obj)
+ return
+
+ if not self.Is_same_dimension(obj) :
+ raise 'Coincidence impossible : objects are not of the same dimension'
+
+ correspond = self.Geometric_coincide(obj)
+ if not correspond :
+ raise 'The objects are not geometrically coincident'
+
+ self.Deep_coincide(obj,correspond)
+ self.relation = obj
+
+
+ def Duplicate(self) :
+
+ return copy.deepcopy(self) # special deepcopy for the Mesh_Descriptor
+
+
+ def Gmsh(self,mesh) :
+
+ if self.num : return # already meshed object
+ self.mesh = mesh # Storing the mesh
+ self.Intermediate_meshing(mesh) # creation of the intermediate objects
+ num = mesh.num+1 # New object number
+ self.Object_meshing(num) # object meshing (with current number num)
+ mesh.num = num # Updating the current gmsh pointer
+ self.num = num # Storing the current object number
+
+
+ def Gmsh_send(self, line_code) :
+
+ self.mesh.command.append(line_code)
+
+
+ def Intermediate_meshing(self,mesh) :
+ pass
+
+
+ def Object_meshing(self,num) :
+ raise "Creation of the Gmsh source not implemented"
+
+
+
+
+# -------------------- POINT OBJECTS ---------------------
+
+
+class Point(Geometric) :
+
+ """
+ POINT OBJECT
+
+ Public methods
+ __init__ :
+ in -> coordinates (the 3rd is zero by default)
+
+ Size : set the size of the neighbouring elements
+ in -> size
+
+ Attractor : define the point as an attractor
+ in -> scale_x : size amplification factor in the x-direction
+ in -> scale_y : size amplification factor in the y-direction
+ in -> distance: influence distance for the perturbation
+
+ Translate : translation of the point
+ in -> x,y,z translation vector (default : z=0)
+
+ Attributes
+ coor : coordinates
+ size : neighbouring element size
+ attractor : parameters of the attractor
+ """
+
+ def __init__(self,x,y,z=0) :
+
+ Geometric.__init__(self)
+ self.coor = array([x,y,z], Float)
+ self.attractor = None
+
+
+ def Geometric_coincide(self,obj) :
+
+ global _CARAC
+ prec = _CARAC['prec']
+
+ d = VectorNorm(self.coor - obj.coor)
+ if d < prec*self.md.size :
+ return 1
+ else :
+ return None
+
+
+ def Size(self,h) :
+
+ self.md.size = float(h)
+
+
+ def Attractor(self, scale_x, scale_y, distance) :
+
+ self.attractor = (float(scale_x), float(scale_y), float(distance))
+
+
+ def Translate(self,x,y,z=0) :
+
+ tran = array([x,y,z]).astype(Float)
+ self.coor = self.coor + tran
+
+
+ def Object_meshing(self,num) :
+
+ ch = (
+ 'Point(' + `num` + ') = {'
+ + `self.coor[0]` + ', '
+ + `self.coor[1]` + ', '
+ + `self.coor[2]` + ', '
+ + `self.md.size` + '};'
+ )
+ self.Gmsh_send(ch)
+
+ if self.attractor :
+ ch = (
+ 'Attractor Point{' + `num` + '} = {'
+ + `self.attractor[0]`+','
+ + `self.attractor[1]`+','
+ + `1./self.attractor[2]` + '};'
+ )
+ self.Gmsh_send(ch)
+
+
+
+# -------------------- LINE OBJECTS ----------------------
+
+
+class Line(Geometric) :
+
+ """
+ LINE OBJECT
+
+
+ Public methods
+
+ Attractor : define the point as an attractor
+ in -> scale_x : size amplification factor in the x-direction
+ in -> scale_y : size amplification factor in the y-direction
+ in -> distance: influence distance for the perturbation
+
+ """
+
+
+ def __init__(self,*points) :
+
+ Geometric.__init__(self)
+
+ if len(points) <=1 :
+ raise "There should be at least two points"
+
+ for point in points :
+ if not point.Is_point() :
+ raise "Arguments should be points"
+
+ self.points = list(points)
+ self.attractor = None
+
+
+ def Geometric_coincide(self,obj) :
+
+ nb_points = len(self.points)
+ if nb_points <> len(obj.points) :
+ raise 'To coincide, lines should have the same number of points'
+
+ # same order of points
+ info = range(nb_points)
+ for i in range(nb_points) :
+ p1 = self.points[i]
+ p2 = obj.points[info[i]]
+ if not p1.Geometric_coincide(p2) :
+ break
+ else :
+ return info
+
+ # reverse order of points
+ info.reverse()
+ for i in range(nb_points) :
+ p1 = self.points[i]
+ p2 = obj.points[info[i]]
+ if not p1.Geometric_coincide(p2) :
+ break
+ else :
+ return info
+
+ return None
+
+
+ def Deep_coincide(self,obj,info) :
+
+ for i in range(len(info)) :
+ p1 = self.points[i]
+ p2 = obj.points[info[i]]
+ p1.Coincide(p2)
+
+
+
+ def Translate(self,x,y,z=0) :
+
+ for point in self.points :
+ point.Translate(x,y,z)
+
+
+ def Transfinite(self,number,progression = 1) :
+
+ self.md.number = int(number)
+ self.md.progression = float(progression)
+
+
+ def Attractor(self,scale_x, scale_y, distance) :
+
+ self.attractor = (float(scale_x), float(scale_y), float(distance))
+
+
+ def __rmul__(self,base) :
+
+ if len(self.points) > 2 :
+ raise "Support (right argument) should be a straight line"
+
+ if self.points[0] in base.points :
+ supp_orig = 0
+ supp_extr = 1
+ elif self.points[1] in base.points :
+ supp_orig = 1
+ supp_extr = 0
+ else :
+ raise "No common point"
+
+ if self.points[supp_orig] == base.points[0] :
+ base_orig = 0
+ base_extr = -1
+ else :
+ base_orig = -1
+ base_extr = 0
+
+ # Translation vector
+ ce = self.points[supp_extr].coor
+ co = self.points[supp_orig].coor
+ tran = ce-co
+
+ # Definition of the edge opposite to the base
+ opp_base = base.Duplicate()
+ opp_base.Translate(tran[0],tran[1],tran[2])
+ opp_base.points[base_orig] = self.points[supp_extr]
+
+ # Definition of the edge opposite to the support
+ opp_supp = self.Duplicate()
+ opp_supp.points[0] = base.points[base_extr]
+ opp_supp.points[1] = opp_base.points[base_extr]
+
+ surf = Surface(base,self,opp_base,opp_supp)
+
+ if len(base.points) > 2 : surf.Ruled()
+
+ return surf
+
+
+ def Intermediate_meshing(self,mesh) :
+
+ for point in self.points :
+ point.Gmsh(mesh)
+
+
+ def Object_meshing(self,num) :
+
+ ch = self.type + '(' + `num` + ') = {'
+ for point in self.points :
+ ch = ch + `point.num` + ','
+ ch = ch[:-1] + '};'
+ self.Gmsh_send(ch)
+
+ if self.md.transfinite :
+ ch = (
+ 'Transfinite Line{' + `num` + '} = ' +
+ `self.md.number+1` +
+ ' Using Progression ' + `self.md.progression` + ';'
+ )
+ self.Gmsh_send(ch)
+
+ if self.attractor :
+ ch = (
+ 'Attractor Line{' + `num` + '} = {'
+ + `self.attractor[0]`+','
+ + `self.attractor[1]`+','
+ + `1./self.attractor[2]` + '};'
+ )
+ self.Gmsh_send(ch)
+
+
+
+class Circle(Line) : pass
+
+# The class inherits everything from Line but its name (to tell Gmsh
+# that it is a circle).
+
+
+
+def Curve(l_x,l_y,l_z=None) :
+
+ if not l_z :
+ l_z = [0.] * len(l_x)
+
+ l_P = []
+ for x,y,z in map(None,l_x,l_y,l_z) :
+ l_P.append(Point(x,y,z))
+
+ line = apply(Line,l_P)
+ return line
+
+
+
+
+
+# -------------------- SURFACE OBJECTS ---------------------
+
+
+class Surface(Geometric) :
+
+ """
+ SURFACE OBJECT (inherit from the Geometric class)
+
+ Public methods
+ __init__ :
+ in -> lines : external bounday of the surface (lines should be connected)
+
+ Holes : set the internal holes (surfaces)
+ in -> holes : list of holes
+
+ Boundary : checks that the boundary is a closed loop and returns the orientation of the edges
+
+ Summit : returns the summit list sorted according to the orientation (see Boundary method)
+
+ Ruled : declare the surface is a ruled one
+
+ Translate : translate the surface
+ in -> x,y,z translation vector (default : z=0)
+
+ Recombine : recombine the surface (try to mesh with quadrangles instead of triangles)
+
+ Transfinite : Declare the mesh to be transfinite
+
+
+ Attributes
+ lines : list of external boundary lines
+ holes : list of internal holes (surfaces)
+ ruled : indicates (false or true) if the surface is a ruled surface
+ loops : list of boundary (external and internal) loops (computed when meshing)
+ """
+
+
+ def __init__(self,*lines) :
+
+ Geometric.__init__(self)
+ self.lines = list(lines)
+ self.holes = []
+ self.ruled = 0
+
+ # Check Assumptions
+ for line in lines :
+ if not line.Is_line() :
+ raise "Arguments should be lines : " + repr(line)
+ if lines == 0 : raise "There should be at least one line"
+ self.Boundary()
+
+
+ def Boundary(self) :
+
+ # checking the boundary is a loop
+ orie = []
+ tmp = list(self.lines) + [self.lines[0]]
+ for i in range(len(self.lines)) :
+ lb = tmp[i]
+ la = tmp[i+1]
+ if lb.points[-1] in [la.points[0], la.points[-1]] :
+ orie.append(1)
+ elif lb.points[0] in [la.points[0], la.points[-1]] :
+ orie.append(-1)
+ else :
+ raise "This is not a loop"
+
+ # checking the boundary is closed
+ if orie[0] == 1 : pi = self.lines[0].points[0]
+ if orie[0] == -1 : pi = self.lines[0].points[-1]
+ if orie[-1] == 1 : pf = self.lines[-1].points[-1]
+ if orie[-1] == -1 : pf = self.lines[-1].points[0]
+ if pi <> pf : raise "The loop is not closed"
+
+ return orie
+
+
+ def Summit(self) :
+
+ summits = []
+ for line, orie in map(None,self.lines,self.Boundary()) :
+ if orie == 1 :
+ summits.append(line.points[0])
+ else :
+ summits.append(line.points[-1])
+ return summits
+
+
+ def Holes(self,*holes) :
+
+ for hole in holes :
+ if not hole.Is_surface() :
+ raise "Holes should be surfaces"
+ self.holes = list(holes)
+
+
+ def Geometric_coincide(self,obj) :
+
+ """
+ return (line_order, hole_order) :
+ line_order : list of the coupled lines ith line of self with line_order[i]th line of obj
+ hole_order : same as line_order but with the internal holes
+ """
+
+ if len(self.lines) <> len(obj.lines) :
+ raise 'To coincide, surfaces should have the same number of border lines'
+
+ if len(self.holes) <> len(obj.holes) :
+ raise 'To coincide, surfaces should have the same number of internal holes'
+
+ # Coincidence of the surface holes
+ hole_order = []
+ nb_holes = len(self.holes)
+ for hole_1 in self.holes :
+ for i in xrange(nb_holes) :
+ if i in hole_order :
+ continue
+ hole_2 = obj.holes[i]
+ if hole_1.Geometric_coincide(hole_2) :
+ hole_order.append(i)
+ break
+ else :
+ return None
+
+ # Coincidence of the external boundary lines
+ line_order = []
+ nb_lines = len(self.lines)
+ for line_1 in self.lines :
+ for i in xrange(nb_lines) :
+ if i in line_order :
+ continue
+ line_2 = obj.lines[i]
+ if line_1.Geometric_coincide(line_2) :
+ line_order.append(i)
+ break
+ else :
+ return None
+
+ return (line_order, hole_order)
+
+
+ def Deep_coincide(self,obj,info) :
+
+ line_order = info[0]
+ hole_order = info[1]
+
+ for i,j in map(None,xrange(len(line_order)),line_order) :
+ l1 = self.lines[i]
+ l2 = obj.lines[j]
+ l1.Coincide(l2)
+
+ for i,j in map(None,xrange(len(hole_order)),hole_order) :
+ h1 = self.holes[i]
+ h2 = obj.holes[j]
+ h1.Coincide(h2)
+
+
+
+ def Ruled(self) :
+
+ self.ruled = 1
+
+ if len(self.lines) not in [3,4] :
+ raise "Ruled surfaces require 3 or 4 edges"
+
+ if self.holes :
+ raise "Holes are forbidden for ruled surfaces"
+
+
+ def Translate(self,x,y,z=0) :
+
+ l_points = []
+ for surf in [self] + self.holes :
+ for line in surf.lines :
+ for point in line.points :
+ if point not in l_points : l_points.append(point)
+
+ for point in l_points :
+ point.Translate(x,y,z)
+
+
+ def Recombine(self,val=1) :
+
+ self.md.recombine = val
+
+
+ def Transfinite(self) :
+
+ self.Ruled()
+
+ if len(self.lines) == 4 :
+ self.Recombine()
+
+ self.md.transfinite = 1
+
+ for line in self.lines :
+ if not line.md.transfinite :
+ raise "Transfinite surfaces require transfinite edges"
+
+ if (
+ self.lines[0].md.number <> self.lines[2].md.number or
+ self.lines[1].md.number <> self.lines[3].md.number
+ ) :
+ raise "Coupled edges should have the same number of elements"
+
+
+
+ def Intermediate_meshing(self,mesh) :
+
+ self.loops = []
+ for surf in [self]+self.holes :
+ loop = LineLoop(surf)
+ self.loops.append(loop)
+ loop.Gmsh(mesh)
+
+
+ def Object_meshing(self,num) :
+
+ # Creation of the surface
+ if self.ruled :
+ ch = 'Ruled Surface(' + `num` + ') = {'
+ else :
+ ch = 'Plane Surface(' + `num` + ') = {'
+ for loop in self.loops :
+ ch = ch + `loop.num` + ','
+ ch = ch[:-1] + '};'
+ self.Gmsh_send(ch)
+
+ # Declaration of transfinite surface
+ if self.md.transfinite :
+ ch = 'Transfinite Surface {' + `num` + '} = {'
+ for summit in self.Summit() :
+ ch = ch + `summit.num` + ','
+ ch = ch[:-1] + '};'
+ self.Gmsh_send(ch)
+
+ # Recombine elements if requested
+ if self.md.recombine :
+ self.Gmsh_send('Recombine Surface {' + `num` + '} ;')
+
+
+
+class LineLoop(Geometric) : # Used only during the meshing phase
+
+
+ def __init__(self,surface) :
+
+ Geometric.__init__(self)
+ self.surface = surface
+
+
+ def Intermediate_meshing(self,mesh) :
+
+ for line in self.surface.lines :
+ line.Gmsh(mesh)
+
+
+ def Object_meshing(self,num) :
+
+ ch = 'Line Loop(' + `num` + ') = {'
+ for line,orie in map(None,self.surface.lines,self.surface.Boundary()) :
+ ch = ch + `orie*line.num` + ','
+ ch = ch[:-1] + '};'
+ self.Gmsh_send(ch)
+
+
+
+
+class Volume(Geometric) :
+
+
+ """
+ VOLUME OBJECT (inherit from the Geometric class)
+
+ Public methods
+ __init__ :
+ in -> surfaces : external bounday of the volume (surfaces should be connected)
+
+ Edge : returns the list of edges
+
+ Holes : set the internal holes (surfaces)
+ in -> holes : list of holes
+
+ Transfinite : Declare the mesh to be transfinite (force the surfaces to be transfinite too)
+
+ Translate : translate the surface
+ in -> x,y,z translation vector (default : z=0)
+
+
+ Attributes
+ surfaces : list of external boundary surfaces
+ holes : list of internal holes (volumes)
+ loops : list of boundary (external and internal) loops (computed when meshing)
+
+
+ Private methods :
+
+ Boundary : checks that the boundary is a closed loop and returns the orientation of the edges
+
+ """
+
+
+ def __init__(self,*surfaces) :
+
+ Geometric.__init__(self)
+ self.surfaces = list(surfaces)
+ self.holes = []
+ self.ruled = 0
+
+ # Check Assumptions
+ for surface in surfaces :
+ if not surface.Is_surface() :
+ raise "Arguments should be surfaces : " + repr(surface)
+ if len(surfaces) < 2 : raise "There should be at least two surfaces"
+ self.Boundary()
+
+
+ def Boundary(self) :
+
+ edges = []
+ for surface in self.surfaces :
+ edges = edges + surface.lines
+
+ # each edge has to appear twice in the list of edges
+ for edge in edges :
+ if edges.count(edge) <> 2 :
+ raise "The surface loop is not closed : each edge should appear twice"
+
+
+ def Edge(self) :
+
+ edges = []
+ for surface in self.surfaces :
+ for line in surface.lines :
+ if line not in edges : edges.append(line)
+ return edges
+
+
+ def Holes(self,*holes) :
+
+ for hole in holes :
+ if not hole.Is_volume() :
+ raise "Holes should be volumes"
+ self.holes = list(holes)
+
+
+ def Geometric_coincide(self,obj) :
+
+ """
+ return (surface_order, hole_order) :
+ surface_order : list of the coupled surfaces ith surface of self with surface_order[i]th surface of obj
+ hole_order : same as surface_order but with the internal holes
+ """
+
+ if len(self.surfaces) <> len(obj.surfaces) :
+ raise 'To coincide, volumes should have the same number of border surfaces'
+
+ if len(self.holes) <> len(obj.holes) :
+ raise 'To coincide, volumes should have the same number of internal holes'
+
+ # Coincidence of the surface holes
+ hole_order = []
+ nb_holes = len(self.holes)
+ for hole_1 in self.holes :
+ for i in xrange(nb_holes) :
+ if i in hole_order :
+ continue
+ hole_2 = obj.holes[i]
+ if hole_1.Geometric_coincide(hole_2) :
+ hole_order.append(i)
+ break
+ else :
+ return None
+
+ # Coincidence of the external boundary lines
+ surface_order = []
+ nb_surfaces = len(self.surfaces)
+ for surface_1 in self.surfaces :
+ for i in xrange(nb_surfaces) :
+ if i in surface_order :
+ continue
+ line_2 = obj.surfaces[i]
+ if surface_1.Geometric_coincide(surface_2) :
+ surface_order.append(i)
+ break
+ else :
+ return None
+
+ return (surface_order, hole_order)
+
+
+ def Deep_coincide(self,obj,info) :
+
+ surface_order = info[0]
+ hole_order = info[1]
+
+ for i,j in map(None,xrange(len(surface_order)),surface_order) :
+ l1 = self.surfaces[i]
+ l2 = obj.surfaces[j]
+ l1.Coincide(l2)
+
+ for i,j in map(None,xrange(len(hole_order)),hole_order) :
+ h1 = self.holes[i]
+ h2 = obj.holes[j]
+ h1.Coincide(h2)
+
+
+ def Transfinite(self) :
+
+ if len(self.surfaces) == 5 :
+ raise "Not implemented"
+
+ if len(self.surfaces) not in [5,6] :
+ raise "Transfinite volumes require 5 or 6 faces"
+
+ if self.holes :
+ raise "Holes are forbidden for transfinite volumes"
+
+ self.md.transfinite = 1
+
+ for surface in self.surfaces :
+ if not surface.md.transfinite :
+ surface.Transfinite() # attention : ce n'est pas vrai dans le cas des prismes
+
+# ATTENTION : ICI, IL FAUDRAIT VERIFIER QUE LES SURFACES PEUVENT ETRE MISES EN VIS A VIS
+
+
+
+ def Translate(self,x,y,z=0) :
+
+ l_points = []
+ for volu in [self] + self.holes :
+ for surf in volu.surfaces :
+ for line in surf.lines :
+ for point in line.points :
+ if point not in l_points : l_points.append(point)
+
+ for point in l_points :
+ point.Translate(x,y,z)
+
+
+
+ def Intermediate_meshing(self,mesh) :
+
+ self.loops = []
+ for volume in [self]+self.holes :
+ loop = SurfaceLoop(volume)
+ self.loops.append(loop)
+ loop.Gmsh(mesh)
+
+
+ def Object_meshing(self,num) :
+
+ # Creation of the volume
+ ch = 'Volume(' + `num` + ') = {'
+ for loop in self.loops :
+ ch = ch + `loop.num` + ','
+ ch = ch[:-1] + '};'
+ self.Gmsh_send(ch)
+
+ # Declaration of transfinite surface
+ if self.md.transfinite :
+
+ bottom_summits = self.surfaces[0].Summit()
+ edges = self.Edge()
+ top_summits = []
+ for summit in bottom_summits :
+ for edge in edges :
+ if summit == edge.points[0] and edge.points[-1] not in bottom_summits :
+ top_summits.append(edge.points[-1])
+ break
+ elif summit == edge.points[-1] and edge.points[0] not in bottom_summits :
+ top_summits.append(edge.points[0])
+ break
+
+ ch = 'Transfinite Volume {' + `num` + '} = {'
+ for summit in bottom_summits + top_summits :
+ ch = ch + `summit.num` + ','
+ ch = ch[:-1] + '};'
+ self.Gmsh_send(ch)
+
+
+
+class SurfaceLoop(Geometric) : # Used only during the meshing phase
+
+
+ def __init__(self,volume) :
+
+ Geometric.__init__(self)
+ self.volume = volume
+
+
+ def Intermediate_meshing(self,mesh) :
+
+ for surface in self.volume.surfaces :
+ surface.Gmsh(mesh)
+
+
+ def Object_meshing(self,num) :
+
+ ch = 'Surface Loop(' + `num` + ') = {'
+ for surface in self.volume.surfaces :
+ ch = ch + `surface.num` + ','
+ ch = ch[:-1] + '};'
+ self.Gmsh_send(ch)
+
+
+
+
+# ------------------- GEOMETRICAL TRANSFORMATION --------------
+
+def VectorProduct(u,v) :
+
+ return array([u[1]*v[2]-u[2]*v[1],u[2]*v[0]-u[0]*v[2],u[0]*v[1]-u[1]*v[0]])
+
+
+def VectorNorm(u) :
+
+ return sqrt(dot(u,u))
+
+
+class Rotation :
+
+ def __init__(self,A,C,B) :
+
+ self.C = c
+ self.a = A-C
+ n = VectorProduct(self.a,B-C)
+ self.n = n / VectorNorm(n)
+
+
+ def Proj(self,M) :
+
+ lbd = dot(M-self.C,self.n)
+ H = self.C + lbd*self.n
+ return H
+
+
+def Scaling_P2(p,t) : return (1.-p)*t*t+p*t
+
+def Scaling_P3(p,t) :
+ q = 1./p
+ a = p+q-2
+ b = 3-2*p-q
+ return a*t**3 + b*t*t +p*t
+
+
+# -------------------- MESHING OPERATIONS ---------------------
+
+class Mesh_Descriptor :
+
+ """
+ Attributes
+ relation Another mesh descriptor provides the mesh parameters
+ parameters dictionnary of the mesh parameters
+ size Point size
+ transfinite Transfinite mesh (0 or 1)
+ number Number of elements along a line (transfinite)
+ progression Progression of element size (transfinite)
+ recombine Recombine mesh or not
+
+ Specific access :
+ md.parameter_name = xxx -> the relation is destroyed (set to None)
+ xxx = md.parameter_name -> if there is a relation, the effective
+ parameter is looked for recursively
+
+ Deep copying : a relation is set to the model instead of a true copy
+ """
+
+ List_Attr = ['size','transfinite','number','progression','recombine']
+
+
+ def __init__(self) :
+
+ self.relation = None
+ self.parameters = {
+ 'size' : 1. , # Point size
+ 'transfinite': 0 , # Transfinite mesh (0 or 1)
+ 'recombine' : 0 # Recombine mesh or not
+ }
+
+
+ def __setattr__(self, attr, value) :
+
+ if attr in Mesh_Descriptor.List_Attr :
+ self.relation = None
+ self.parameters[attr] = value
+
+ if attr == 'number' :
+ self.transfinite = 1
+
+ else :
+ self.__dict__[attr] = value
+
+
+ def __getattr__(self,attr) :
+
+ if self.relation :
+ return (getattr(self.relation,attr))
+ else :
+ if attr in self.parameters.keys() :
+ return self.parameters[attr]
+ else :
+ raise AttributeError
+
+
+ def __deepcopy__(self,visit) :
+
+ md = copy.copy(self)
+ md.parameters = copy.copy(self.parameters)
+ md.relation = self
+ return md
+
+
+
+class Mesh :
+
+ """
+
+ """
+
+
+ def __init__(self, algo = 2, order = 1, gmsh='gmsh') :
+
+ self.num_ph = 0
+ self.num = 0
+ self.order = order
+ self.command = ['Mesh.Algorithm = ' + repr(algo) + ' ;']
+ self.command += ['Mesh.ElementOrder = ' + repr(order) + ' ;']
+ self.physicals = {}
+ self.gmsh = gmsh
+
+
+ def Physical(self, name, *l_lobj) :
+
+ # Checking the name
+ if type(name) <> type(' ') :
+ raise 'First argument should be the name of the physical'
+ if name in self.physicals.keys() :
+ raise 'Physical '+name+' already exists'
+
+ # treating the case of list of lists parameters
+ l_obj = []
+ for l in l_lobj :
+ if type(l) == type([]) :
+ l_obj = l_obj + l
+ else :
+ l_obj.append(l)
+
+ # Checking all objects are geometric
+ for obj in l_obj :
+ if not Is_Geometric(obj) :
+ raise "Non geometrical object : " + repr(obj) + " Physical = " + name
+
+ cl = l_obj[0].Base_class()
+ # Checking all objects are of the same dimension
+ # ref_dim = l_obj[0]
+ # for obj in l_obj[1:] :
+ # if not ref_dim.Is_same_dimension(obj) :
+ # raise "All objects are not of the same dimension : " + repr(obj)
+
+ # Creation of the objects if necessary
+ for obj in l_obj :
+ obj.Gmsh(self)
+
+ # Creation of the physical
+ self.num_ph= self.num_ph + 1
+ ch = name + '=' + `self.num_ph` + ';'
+ self.command.append(ch)
+ ch = 'Physical ' + cl + '(' + name + ') = {'
+ for obj in l_obj :
+ ch = ch + `obj.num` + ','
+ ch = ch[:-1] + '};'
+ self.command.append(ch)
+
+ # Name of the physical
+ name_gmsh = 'GM'+`self.num_ph`
+ self.physicals[name] = name_gmsh
+
+
+ def Save(self, file = 'fort.geo') :
+
+ if os.path.isfile(file) :
+ os.remove(file)
+
+ f = open(file,'w')
+ f.write(string.joinfields(self.command,'\n'))
+ f.close()
+
+
+ def View(self) :
+
+ self.Save('fort.geo')
+# os.system('gmsh fort.geo')
+ os.system(self.gmsh + ' fort.geo')
+ os.remove('fort.geo')
+
+
+ def Create(self, file = 'fort.19') :
+
+ self.Save()
+# os.system('gmsh -3 fort.geo')
+ os.system(self.gmsh + ' -3 fort.geo')
+ os.rename('fort.msh',file)
+
+
+ def Name(self, MA, CREA_GROUP_NO) :
+
+ l_gma = []
+ l_mcf = []
+ for gma in self.physicals.keys() :
+ l_gma.append(self.physicals[gma])
+ l_mcf.append(_F(GROUP_MA = self.physicals[gma],NOM=gma))
+
+ DEFI_GROUP(reuse = MA,
+ MAILLAGE = MA,
+ CREA_GROUP_MA = tuple(l_mcf),
+ )
+
+ SMESH_02 = CREA_MAILLAGE(
+ MAILLAGE = MA,
+ DETR_GROUP_MA = _F(GROUP_MA = tuple(l_gma)),
+ )
+
+ DETRUIRE(CONCEPT = _F(NOM = MA), INFO=1)
+
+ if CREA_GROUP_NO == 'OUI' :
+ DEFI_GROUP(reuse = SMESH_02,
+ MAILLAGE = SMESH_02,
+ CREA_GROUP_NO = _F(TOUT_GROUP_MA = 'OUI'),
+ )
+
+ else :
+# Traitement des GROUP_NO qui sont des points
+ info_gno = SMESH_02.LIST_GROUP_NO()
+ l_gno = []
+ for gno in info_gno :
+ if gno[1] == 1 : l_gno.append(gno[0])
+
+ l_gma = []
+ for gma in self.physicals.keys() :
+ nom_gmsh = self.physicals[gma]
+ if nom_gmsh in l_gno :
+ l_gma.append(gma)
+
+ if l_gma :
+ DEFI_GROUP(reuse = SMESH_02,
+ MAILLAGE = SMESH_02,
+ CREA_GROUP_NO = _F(GROUP_MA = tuple(l_gma)),
+ )
+
+ return SMESH_02
+
+
+
+ def LIRE_GMSH(self,
+ UNITE_GMSH = 19,
+ UNITE_MAILLAGE = 20,
+ MODI_QUAD = 'NON',
+ CREA_GROUP_NO = 'OUI'
+ ) :
+
+ """
+ Lecture du maillage (format Aster) a partir de sa definition
+ (format sup_gmsh)
+ UNITE_GMSH = Numero d'unite logique pour le fichier msh
+ UNITE_MAILLAGE = Numero d'unite logique pour le fichier mail
+ MODI_QUAD = 'OUI' si line->quad, 'NON' sinon
+ CREA_GROUP_NO = 'OUI' si on cree les group_no, 'NON' sinon
+ """
+
+ nom_gmsh = 'fort.' + repr(UNITE_GMSH)
+ self.Create(nom_gmsh)
+
+ PRE_GMSH(UNITE_GMSH=UNITE_GMSH, UNITE_MAILLAGE=UNITE_MAILLAGE)
+
+ SMESH_00 = LIRE_MAILLAGE(UNITE = UNITE_MAILLAGE)
+ DEFI_FICHIER(ACTION='LIBERER',UNITE = UNITE_GMSH)
+ DEFI_FICHIER(ACTION='LIBERER',UNITE = UNITE_MAILLAGE)
+
+ if MODI_QUAD == 'OUI' and self.order == 2 :
+ raise 'The finite elements are already of second order'
+
+ if MODI_QUAD == 'OUI' and self.order <> 2 :
+ SMESH_01 = CREA_MAILLAGE(
+ MAILLAGE = SMESH_00,
+ LINE_QUAD = _F(TOUT = 'OUI')
+ )
+ DETRUIRE(CONCEPT=_F(NOM=SMESH_00), INFO=1)
+ SMESH_00 = SMESH_01
+
+ SMESH_00 = self.Name(SMESH_00,CREA_GROUP_NO)
+
+ return SMESH_00
--- /dev/null
+#@ MODIF t_fonction Utilitai DATE 31/05/2005 AUTEUR DURAND C.DURAND
+# -*- coding: iso-8859-1 -*-
+# CONFIGURATION MANAGEMENT OF EDF VERSION
+# ======================================================================
+# COPYRIGHT (C) 1991 - 2005 EDF R&D WWW.CODE-ASTER.ORG
+# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
+# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
+# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
+# (AT YOUR OPTION) ANY LATER VERSION.
+#
+# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
+# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
+# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
+# GENERAL PUBLIC LICENSE FOR MORE DETAILS.
+#
+# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
+# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
+# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
+# ======================================================================
+from Numeric import *
+import copy
+import types
+
+def interp(typ_i,val,x1,x2,y1,y2) :
+ if typ_i==['LIN','LIN']: return y1+(y2-y1)*(val-x1)/(x2-x1)
+ if typ_i==['LIN','LOG']: return exp(log(y1)+(val-x1)*(log(y2)-log(y1))/(x2-x1))
+ if typ_i==['LOG','LOG']: return exp(log(y1)+(log(val)-log(x1))*(log(y2)-log(y1))/(log(x2)-log(x1)))
+ if typ_i==['LOG','LIN']: return y1+(log(val)-log(x1))*(y2-y1)/(log(x2)-log(x1))
+ if typ_i[0]=='NON' :
+ if val==x1 : return y1
+ elif val==x2 : return y2
+ else : raise StandardError, 'fonction : interpolation NON'
+def is_ordo(liste) :
+ listb=dict([(i,0) for i in liste]).keys()
+ listb.sort()
+ return liste==listb
+
+class t_fonction :
+ ### Classe pour fonctions réelles, équivalent au type aster = fonction_sdaster
+ def __init__(self,vale_x,vale_y,para) :
+ # création d'un objet fonction
+ # vale_x et vale_y sont des listes de réels de meme longueur
+ # para est un dictionnaire contenant les entrées PROL_DROITE, PROL_GAUCHE et INTERPOL (cf sd ASTER)
+ pk=para.keys()
+ pk.sort()
+ if pk!=['INTERPOL','NOM_PARA','NOM_RESU','PROL_DROITE','PROL_GAUCHE'] :
+ raise StandardError, 'fonction : parametres incorrects'
+ if para['INTERPOL'] not in [['NON','NON'],['LIN','LIN'],['LIN','LOG'],['LOG','LOG'],['LOG','LIN'],] :
+ raise StandardError, 'fonction : parametre INTERPOL incorrect'
+ if para['PROL_DROITE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
+ raise StandardError, 'fonction : parametre PROL_DROITE incorrect'
+ if para['PROL_GAUCHE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
+ raise StandardError, 'fonction : parametre PROL_GAUCHE incorrect'
+ self.vale_x = array(vale_x)
+ self.vale_y = array(vale_y)
+ self.para = para
+ if len(self.vale_x)!=len(self.vale_y) :
+ raise StandardError, 'fonction : longueur abscisse <> longueur ordonnées'
+ if not is_ordo(self.vale_x) :
+ raise StandardError, 'fonction : abscisses non strictement croissantes'
+
+ def __add__(self,other) :
+ # addition avec une autre fonction ou un nombre, par surcharge de l'opérateur +
+ if isinstance(other,t_fonction):
+ para=copy.copy(self.para)
+ vale_x,para['PROL_GAUCHE'],para['PROL_DROITE']=self.homo_support(other)
+ fff=self.evalfonc(vale_x)
+ ggg=other.evalfonc(vale_x)
+ if isinstance(self,t_fonction_c): return t_fonction_c(vale_x,fff.vale_y+ggg.vale_y,para)
+ else : return t_fonction(vale_x,fff.vale_y+ggg.vale_y,para)
+ elif type(other) in [types.FloatType,types.IntType,types.ComplexType] :
+ if isinstance(self,t_fonction_c): return t_fonction_c(self.vale_x,self.vale_y+other,self.para)
+ else : return t_fonction(self.vale_x,self.vale_y+other,self.para)
+ else: raise StandardError, 'fonctions : erreur de type dans __add__'
+
+ def __mul__(self,other) :
+ # multiplication avec une autre fonction ou un nombre, par surcharge de l'opérateur *
+ if isinstance(other,t_fonction):
+ para=copy.copy(self.para)
+ vale_x,para['PROL_GAUCHE'],para['PROL_DROITE']=self.homo_support(other)
+ fff=self.evalfonc(vale_x)
+ ggg=other.evalfonc(vale_x)
+ if isinstance(self,t_fonction_c): return t_fonction_c(vale_x,fff.vale_y*ggg.vale_y,para)
+ else : return t_fonction(vale_x,fff.vale_y*ggg.vale_y,para)
+ elif type(other) in [types.FloatType,types.IntType] :
+ return t_fonction(self.vale_x,self.vale_y*other,self.para)
+ elif type(other) ==types.ComplexType :
+ return t_fonction_c(self.vale_x,self.vale_y*other,self.para)
+ else: raise StandardError, 'fonctions : erreur de type dans __mul__'
+
+ def __repr__(self) :
+ # affichage de la fonction en double colonne
+ texte=[]
+ for i in range(len(self.vale_x)) :
+ texte.append('%f %f' % (self.vale_x[i],self.vale_y[i]))
+ return '\n'.join(texte)
+
+ def __getitem__(self,other) :
+ # composition de deux fonction F[G]=FoG=F(G(x))
+ para=copy.copy(self.para)
+ if other.para['NOM_RESU']!=self.para['NOM_PARA'] :
+ raise StandardError,'''composition de fonctions : NOM_RESU1 et NOM_PARA2 incohérents '''
+ para['NOM_PARA']==other.para['NOM_PARA']
+ return t_fonction(other.vale_x,map(self,other.vale_y),para)
+
+ def __call__(self,val,tol=1.e-6):
+ # méthode pour évaluer f(x)
+ # tolérance, par défaut 1.e-6 en relatif sur la longueur de l'intervalle
+ # adjacent, pour capter les erreurs d'arrondi en cas de prolongement exclu
+ i=searchsorted(self.vale_x,val)
+ n=len(self.vale_x)
+ if i==0 :
+ if self.para['PROL_GAUCHE']=='EXCLU' :
+ eps_g=(val-self.vale_x[0] )/(self.vale_x[1] -self.vale_x[0])
+ if abs(eps_g)<=tol : return self.vale_y[0]
+ else : raise StandardError, 'fonction évaluée hors du domaine de définition'
+ else :
+ if self.para['PROL_GAUCHE']=='CONSTANT' : return self.vale_y[0]
+ if self.para['PROL_GAUCHE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val,self.vale_x[0],
+ self.vale_x[1],
+ self.vale_y[0],
+ self.vale_y[1])
+ elif i==n :
+ if self.para['PROL_DROITE']=='EXCLU' :
+ eps_d=(val-self.vale_x[-1])/(self.vale_x[-1]-self.vale_x[-2])
+ if abs(eps_d)<=tol : return self.vale_y[-1]
+ else : raise StandardError, 'fonction évaluée hors du domaine de définition'
+ else :
+ if self.para['PROL_DROITE']=='CONSTANT' : return self.vale_y[-1]
+ if self.para['PROL_DROITE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val,self.vale_x[-1],
+ self.vale_x[-2],
+ self.vale_y[-1],
+ self.vale_y[-2])
+ else :
+ return interp(self.para['INTERPOL'],val,self.vale_x[i-1],
+ self.vale_x[i],
+ self.vale_y[i-1],
+ self.vale_y[i])
+
+ def homo_support(self,other) :
+ # renvoie le support d'abscisses homogénéisé entre self et other
+ # i.e. si prolongement exclu, on retient plus grand min ou plus petit max, selon
+ # si prolongement autorisé, on conserve les abscisses d'une fonction, extrapolantes
+ # sur l'autre.
+ # Pour les points intermédiaires : union et tri des valeurs des vale_x réunis.
+ if other.vale_x[0]>self.vale_x[0]:
+ if other.para['PROL_GAUCHE']!='EXCLU' : f_g=self
+ else : f_g=other
+ else :
+ if self.para['PROL_GAUCHE'] !='EXCLU' : f_g=other
+ else : f_g=self
+ val_min =f_g.vale_x[0]
+ prol_gauche=f_g.para['PROL_GAUCHE']
+ if self.vale_x[0]>other.vale_x[0]:
+ if other.para['PROL_DROITE']!='EXCLU' : f_d=self
+ else : f_d=other
+ else :
+ if self.para['PROL_DROITE'] !='EXCLU' : f_d=other
+ else : f_d=self
+ val_max =f_d.vale_x[-1]
+ prol_droite=f_d.para['PROL_DROITE']
+ vale_x=concatenate((self.vale_x,other.vale_x))
+ vale_x=clip(vale_x,val_min,val_max)
+ vale_x=sort(dict([(i,0) for i in vale_x]).keys())
+ return vale_x,prol_gauche,prol_droite
+
+ def cut(self,rinf,rsup,prec,crit='RELATIF') :
+ # renvoie la fonction self dont on a 'coupé' les extrémités en x=rinf et x=rsup
+ # pour la recherche de rinf et rsup dans la liste d'abscisses :
+ # prec=precision crit='absolu' ou 'relatif'
+ para=copy.copy(self.para)
+ para['PROL_GAUCHE']='EXCLU'
+ para['PROL_DROITE']='EXCLU'
+ if crit=='ABSOLU' : rinf_tab=greater(abs(self.vale_x-rinf),prec)
+ elif crit=='RELATIF': rinf_tab=greater(abs(self.vale_x-rinf),prec*rinf)
+ else : raise StandardError, 'fonction : cut : critère absolu ou relatif'
+ if crit=='ABSOLU' : rsup_tab=greater(abs(self.vale_x-rsup),prec)
+ elif crit=='RELATIF': rsup_tab=greater(abs(self.vale_x-rsup),prec*rsup)
+ else : raise StandardError, 'fonction : cut : critère absolu ou relatif'
+ if alltrue(rinf_tab) : i=searchsorted(self.vale_x,rinf)
+ else : i=rinf_tab.tolist().index(0)+1
+ if alltrue(rsup_tab) : j=searchsorted(self.vale_x,rsup)
+ else : j=rsup_tab.tolist().index(0)
+ vale_x=array([rinf,]+self.vale_x.tolist()[i:j]+[rsup,])
+ vale_y=array([self(rinf),]+self.vale_y.tolist()[i:j]+[self(rsup),])
+ return t_fonction(vale_x,vale_y,para)
+
+ def cat(self,other,surcharge) :
+ # renvoie une fonction concaténée avec une autre, avec règles de surcharge
+ para=copy.copy(self.para)
+ if self.para['INTERPOL']!=other.para['INTERPOL'] : raise StandardError, 'concaténation de fonctions à interpolations différentes'
+ if min(self.vale_x)<min(other.vale_x) :
+ f1=self
+ f2=other
+ else :
+ f1=other
+ f2=self
+ para['PROL_GAUCHE']=f1.para['PROL_GAUCHE']
+ if surcharge=='GAUCHE' :
+ i=searchsorted(f2.vale_x,f1.vale_x[-1])
+ if i!=len(f2.vale_x) : para['PROL_DROITE']=f2.para['PROL_DROITE']
+ else : para['PROL_DROITE']=f1.para['PROL_DROITE']
+ vale_x=array(f1.vale_x.tolist()+f2.vale_x[i:].tolist())
+ vale_y=array(f1.vale_y.tolist()+f2.vale_y[i:].tolist())
+ elif surcharge=='DROITE' :
+ i=searchsorted(f1.vale_x,f2.vale_x[0])
+ if i!=len(f1.vale_x) : para['PROL_DROITE']=f2.para['PROL_DROITE']
+ else : para['PROL_DROITE']=f1.para['PROL_DROITE']
+ vale_x=array(f1.vale_x[:i].tolist()+f2.vale_x.tolist())
+ vale_y=array(f1.vale_y[:i].tolist()+f2.vale_y.tolist())
+ return t_fonction(vale_x,vale_y,para)
+
+ def tabul(self) :
+ # mise en forme de la fonction selon un vecteur unique (x1,y1,x2,y2,...)
+ __tab=array([self.vale_x,self.vale_y])
+ return ravel(transpose(__tab)).tolist()
+
+ def extreme(self) :
+ # renvoie un dictionnaire des valeurs d'ordonnées min et max
+ val_min=min(self.vale_y)
+ val_max=max(self.vale_y)
+ vm={}
+ vm['min']=[[self.vale_y[i],self.vale_x[i]] for i in range(len(self.vale_x))\
+ if self.vale_y[i]==val_min]
+ vm['max']=[[self.vale_y[i],self.vale_x[i]] for i in range(len(self.vale_x))\
+ if self.vale_y[i]==val_max]
+ vm['min'].sort()
+ vm['max'].sort()
+ for item in vm['min'] : item.reverse()
+ for item in vm['max'] : item.reverse()
+ return vm
+
+ def trapeze(self,coef) :
+ # renvoie la primitive de la fonction, calculée avec la constante d'intégration 'coef'
+ trapz = zeros(len(self.vale_y),Float)
+ trapz[0] = coef
+ trapz[1:] = (self.vale_y[1:]+self.vale_y[:-1])/2*(self.vale_x[1:]-self.vale_x[:-1])
+ prim_y=cumsum(trapz)
+ para=copy.copy(self.para)
+ para['PROL_GAUCHE']='EXCLU'
+ para['PROL_DROITE']='EXCLU'
+ if para['NOM_RESU'][:4]=='VITE' : para['NOM_RESU']='DEPL'
+ elif para['NOM_RESU'][:4]=='ACCE' : para['NOM_RESU']='VITE'
+ return t_fonction(self.vale_x,prim_y,para)
+
+ def simpson(self,coef) :
+ # renvoie la primitive de la fonction, calculée avec la constante d'intégration 'coef'
+ para=copy.copy(self.para)
+ para['PROL_GAUCHE']='EXCLU'
+ para['PROL_DROITE']='EXCLU'
+ if para['NOM_RESU'][:4]=='VITE' : para['NOM_RESU']='DEPL'
+ elif para['NOM_RESU'][:4]=='ACCE' : para['NOM_RESU']='VITE'
+ fm = self.vale_y[0]
+ fb = self.vale_y[1]
+ h2 = self.vale_x[1] - self.vale_x[0]
+ tabl = [coef,coef +(fb+fm)*h2/2.]
+ prim_y = copy.copy(tabl)
+ iperm = 0
+ ip = (1,0)
+ eps = 1.E-4
+ for i in range(2,len(self.vale_x)) :
+ h1 = h2
+ h2 = self.vale_x[i] - self.vale_x[i-1]
+ bma = h1 + h2
+ fa = fm
+ fm = fb
+ fb = self.vale_y[i]
+ deltah = h2 - h1
+ if h1==0. or h2==0. or abs( deltah / h1 ) <= eps :
+ ct = (1.,4.,1.)
+ else :
+ epsi = deltah / (h1*h2)
+ ct = (1.-epsi*h2,2.+epsi*deltah,1.+epsi*h1)
+ tabl[iperm] = tabl[iperm] + (bma/6.)*(ct[0]*fa+ct[1]*fm+ct[2]*fb)
+ prim_y.append(tabl[iperm])
+ iperm = ip[iperm]
+ return t_fonction(self.vale_x,prim_y,para)
+
+ def derive(self) :
+ # renvoie la dérivée de la fonction
+ pas=self.vale_x[1:]-self.vale_x[:-1]
+ pentes=(self.vale_y[1:]-self.vale_y[:-1])/(self.vale_x[1:]-self.vale_x[:-1])
+ derive=(pentes[1:]*pas[1:]+pentes[:-1]*pas[:-1])/(pas[1:]+pas[:-1])
+ derv_y=[pentes[0]]+derive.tolist()+[pentes[-1]]
+ para=copy.copy(self.para)
+ para['PROL_GAUCHE']='EXCLU'
+ para['PROL_DROITE']='EXCLU'
+ if para['NOM_RESU'][:4]=='DEPL' : para['NOM_RESU']='VITE'
+ elif para['NOM_RESU'][:4]=='VITE' : para['NOM_RESU']='ACCE'
+ return t_fonction(self.vale_x,derv_y,para)
+
+ def inverse(self) :
+ # renvoie l'inverse de la fonction
+ # on intervertit vale_x et vale_y, on swape interpolation
+ para=copy.copy(self.para)
+ para['NOM_RESU']='TOUTRESU'
+ para['NOM_PARA']=self.para['NOM_PARA']
+ para['INTERPOL'].reverse()
+ if para['PROL_GAUCHE']=='CONSTANT' : para['PROL_GAUCHE']='EXCLU'
+ if para['PROL_DROITE']=='CONSTANT' : para['PROL_DROITE']='EXCLU'
+ vale_x=self.vale_y
+ vale_y=self.vale_x
+ if not is_ordo(vale_x) :
+ vale_x=vale_x[::-1]
+ vale_y=vale_y[::-1]
+ return t_fonction(vale_x,vale_y,para)
+
+ def abs(self) :
+ # renvoie la mm fonction avec valeur absolue des ordonnées
+ para=copy.copy(self.para)
+ if para['PROL_GAUCHE']=='LINEAIRE' : para['PROL_GAUCHE']='EXCLU'
+ if para['PROL_DROITE']=='LINEAIRE' : para['PROL_DROITE']='EXCLU'
+ return t_fonction(self.vale_x,absolute(self.vale_y),para)
+
+ def evalfonc(self,liste_val) :
+ # renvoie la mm fonction interpolée aux points définis par la liste 'liste_val'
+ return t_fonction(liste_val,map(self,liste_val),self.para)
+
+ def sup(self,other) :
+ # renvoie l'enveloppe supérieure de self et other
+ para=copy.copy(self.para)
+# commentaire : pour les prolongements et l'interpolation, c'est self
+# qui prime sur other
+ vale_x=self.vale_x.tolist()+other.vale_x.tolist()
+# on ote les abscisses doublons
+ vale_x=dict([(i,0) for i in vale_x]).keys()
+ vale_x.sort()
+ vale_x=array(vale_x)
+ vale_y1=map(self, vale_x)
+ vale_y2=map(other,vale_x)
+ vale_y=map(max,vale_y1,vale_y2)
+ return t_fonction(vale_x,vale_y,para)
+
+ def inf(self,other) :
+ # renvoie l'enveloppe inférieure de self et other
+ para=copy.copy(self.para)
+# commentaire : pour les prolongements et l'interpolation, c'est self
+# qui prime sur other
+ vale_x=self.vale_x.tolist()+other.vale_x.tolist()
+# on ote les abscisses doublons
+ vale_x=dict([(i,0) for i in vale_x]).keys()
+ vale_x.sort()
+ vale_x=array(vale_x)
+ vale_y1=map(self, vale_x)
+ vale_y2=map(other,vale_x)
+ vale_y=map(min,vale_y1,vale_y2)
+ return t_fonction(vale_x,vale_y,para)
+
+ def suppr_tend(self) :
+ # pour les corrections d'accélérogrammes
+ # suppression de la tendance moyenne d'une fonction
+ para=copy.copy(self.para)
+ xy=sum(self.vale_x*self.vale_y)
+ x0=sum(self.vale_x)
+ y0=sum(self.vale_y)
+ xx=sum(self.vale_x*self.vale_x)
+ n=len(self.vale_x)
+ a1 = ( n*xy - x0*y0) / (n*xx - x0*x0)
+ a0 = (xx*x0 - x0*xy) / (n*xx - x0*x0)
+ return t_fonction(self.vale_x,self.vale_y-a1*self.vale_x-a0,self.para)
+
+ def normel2(self) :
+ # norme de la fonction
+ __ex=self*self
+ __ex=__ex.trapeze(0.)
+ return sqrt(__ex.vale_y[-1])
+
+ def fft(self,methode) :
+ # renvoie la transformée de Fourier rapide FFT
+ import FFT
+ para=copy.copy(self.para)
+ para['NOM_PARA']='FREQ'
+ if self.para['NOM_PARA']!='INST' :
+ raise StandardError, 'fonction réelle : FFT : NOM_PARA=INST pour une transformée directe'
+ pas = self.vale_x[1]-self.vale_x[0]
+ for i in range(1,len(self.vale_x)) :
+ ecart = abs(((self.vale_x[i]-self.vale_x[i-1])-pas)/pas)
+ if ecart>1.e-2 :
+ raise StandardError, 'fonction réelle : FFT : la fonction doit etre à pas constant'
+ n=int(log(len(self.vale_x))/log(2))
+ if methode=='TRONCATURE' :
+ vale_y=self.vale_y[:2**n]
+ elif methode=='PROL_ZERO' :
+ vale_y=self.vale_y.tolist()
+ vale_y=vale_y+[0.]*(2**(n+1)-len(self.vale_x))
+ vale_y=array(vale_y)
+ vect=FFT.fft(vale_y)
+ pasfreq =1./(pas*(len(vect)-1))
+ vale_x =[pasfreq*i for i in range(len(vect))]
+ vale_y =vect*pas
+ return t_fonction_c(vale_x,vale_y,para)
+
+class t_fonction_c(t_fonction) :
+ ### Classe pour fonctions complexes, équivalent au type aster = fonction_c
+ def tabul(self) :
+ # mise en forme de la fonction selon un vecteur unique (x1,yr1,yi1,x2,yr2,yr2,...)
+ __tab=array([self.vale_x,self.vale_y.real,self.vale_y.imag])
+ return ravel(transpose(__tab)).tolist()
+
+ def __repr__(self) :
+ # affichage de la fonction en double colonne
+ texte=[]
+ for i in range(len(self.vale_x)) :
+ texte.append('%f %f + %f .j' % (self.vale_x[i],self.vale_y[i].real,self.vale_y[i].imag))
+ return '\n'.join(texte)
+
+ def fft(self,methode,syme) :
+ # renvoie la transformée de Fourier rapide FFT (sens inverse)
+ import FFT
+ para=copy.copy(self.para)
+ para['NOM_PARA']='INST'
+ if self.para['NOM_PARA']!='FREQ' :
+ raise StandardError, 'fonction complexe : FFT : NOM_PARA=FREQ pour une transformée directe'
+ pas = self.vale_x[1]-self.vale_x[0]
+ for i in range(1,len(self.vale_x)) :
+ ecart = abs(((self.vale_x[i]-self.vale_x[i-1])-pas)/pas)
+ if ecart>1.e-3 :
+ raise StandardError, 'fonction complexe : FFT : la fonction doit etre à pas constant'
+ n=int(log(len(self.vale_x))/log(2))
+ if syme=='OUI' and len(self.vale_x)==2**n :
+ vale_fonc=self.vale_y
+ elif syme=='NON' and len(self.vale_x)==2**n :
+ vale_fonc=self.vale_y.tolist()
+ vale_fon1=vale_fonc[:]
+ vale_fon1.reverse()
+ vale_fonc=vale_fonc+vale_fon1
+ vale_fonc=array(vale_fonc)
+ elif syme=='NON' and len(self.vale_x)!=2**n and methode=='PROL_ZERO' :
+ vale_fonc=self.vale_y.tolist()+[complex(0.)]*(2**(n+1)-len(self.vale_x))
+ vale_fon1=vale_fonc[:]
+ vale_fon1.reverse()
+ vale_fonc=vale_fonc+vale_fon1
+ vale_fonc=array(vale_fonc)
+ elif syme=='NON' and len(self.vale_x)!=2**n and methode=='TRONCATURE' :
+ vale_fonc=self.vale_y[:2**n]
+ vale_fonc=vale_fonc.tolist()
+ vale_fon1=vale_fonc[:]
+ vale_fon1.reverse()
+ vale_fonc=vale_fonc+vale_fon1
+ vale_fonc=array(vale_fonc)
+ if syme=='OUI' and len(self.vale_x)!=2**n :
+ raise StandardError, 'fonction complexe : FFT : syme=OUI et nombre de points<>2**n'
+ part1=vale_fonc[:len(vale_fonc)/2+1]
+ part2=vale_fonc[1:len(vale_fonc)/2]
+ part2=conjugate(part2)
+ part2=part2.tolist()
+ part2.reverse()
+ vale_fonc=array(part1.tolist()+part2)
+ vect=FFT.inverse_fft(vale_fonc)
+ vect=vect.real
+ pasfreq =1./(pas*(len(vect)-1))
+ vale_x =[pasfreq*i for i in range(len(vect))]
+ pas2 =(1./self.vale_x[-1])*((len(self.vale_x))/float(len(vect)))
+ vale_y =vect/pas2
+ return t_fonction(vale_x,vale_y,para)
+
+
+class t_nappe :
+ ### Classe pour nappes, équivalent au type aster = nappe_sdaster
+ def __init__(self,vale_para,l_fonc,para) :
+ # création d'un objet nappe
+ # vale_para est la liste de valeur des parametres (mot clé PARA dans DEFI_NAPPE)
+ # para est un dictionnaire contenant les entrées PROL_DROITE, PROL_GAUCHE et INTERPOL (cf sd ASTER)
+ # l_fonc est la liste des fonctions, de cardinal égal à celui de vale_para
+ pk=para.keys()
+ pk.sort()
+ if pk!=['INTERPOL','NOM_PARA','NOM_PARA_FONC','NOM_RESU','PROL_DROITE','PROL_GAUCHE'] :
+ raise StandardError, 'nappe : parametres incorrects'
+ if para['INTERPOL'] not in [['NON','NON'],['LIN','LIN'],
+ ['LIN','LOG'],['LOG','LOG'],['LOG','LIN'],] :
+ raise StandardError, 'nappe : parametre INTERPOL incorrect'
+ if para['PROL_DROITE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
+ raise StandardError, 'nappe : parametre PROL_DROITE incorrect'
+ if para['PROL_GAUCHE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
+ raise StandardError, 'nappe : parametre PROL_GAUCHE incorrect'
+ self.vale_para = array(vale_para)
+ if type(l_fonc) not in (types.ListType,types.TupleType) :
+ raise StandardError, 'nappe : la liste de fonctions fournie n est pas une liste'
+ if len(l_fonc)!=len(vale_para) :
+ raise StandardError, 'nappe : nombre de fonctions différent du nombre de valeurs du paramètre'
+ for f in l_fonc :
+ if not isinstance(f,t_fonction) and not isinstance(f,t_fonction_c) :
+ raise StandardError, 'nappe : les fonctions fournies ne sont pas du bon type'
+ self.l_fonc = l_fonc
+ self.para = para
+
+ def __call__(self,val1,val2):
+ # méthode pour évaluer nappe(val1,val2)
+ i=searchsorted(self.vale_para,val1)
+ n=len(self.vale_para)
+ if i==0 :
+ if val1==self.vale_para[0] : return self.l_fonc[0](val2)
+ if val1 <self.vale_para[0] :
+ if self.para['PROL_GAUCHE']=='EXCLU' : raise StandardError, 'nappe évaluée hors du domaine de définition'
+ if self.para['PROL_GAUCHE']=='CONSTANT' : return self.l_fonc[0](val2)
+ if self.para['PROL_GAUCHE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val1,
+ self.vale_para[0],
+ self.vale_para[1],
+ self.l_fonc[0](val2),
+ self.l_fonc[1](val2))
+ elif i==n :
+ if val1==self.vale_para[-1] : return self.l_fonc[-1](val2)
+ if val1 >self.vale_para[-1] :
+ if self.para['PROL_DROITE']=='EXCLU' : raise StandardError, 'nappe évaluée hors du domaine de définition'
+ if self.para['PROL_DROITE']=='CONSTANT' : return self.l_fonc[-1](val2)
+ if self.para['PROL_DROITE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val1,
+ self.vale_para[-1],
+ self.vale_para[-2],
+ self.l_fonc[-1](val2),
+ self.l_fonc[-2](val2))
+ else :
+ return interp(self.para['INTERPOL'],val1,self.vale_para[i-1],
+ self.vale_para[i],
+ self.l_fonc[i-1](val2),
+ self.l_fonc[i](val2))
+
+ def __add__(self,other) :
+ # addition avec une autre nappe ou un nombre, par surcharge de l'opérateur +
+ l_fonc=[]
+ if isinstance(other,t_nappe):
+ if self.vale_para!=other.vale_para : raise StandardError, 'nappes à valeurs de paramètres différentes'
+ for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]+other.l_fonc[i])
+ elif type(other) in [types.FloatType,types.IntType] :
+ for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]+other)
+ else: raise StandardError, 't_nappe : erreur de type dans __add__'
+ return t_nappe(self.vale_para,l_fonc,self.para)
+
+ def __mul__(self,other) :
+ # multiplication avec une autre fonction ou un nombre, par surcharge de l'opérateur *
+ l_fonc=[]
+ if isinstance(other,t_nappe):
+ if self.vale_para!=other.vale_para : raise StandardError, 'nappes à valeurs de paramètres différentes'
+ for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]*other.l_fonc[i])
+ elif type(other) in [types.FloatType,types.IntType] :
+ for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]*other)
+ else: raise StandardError, 't_nappe : erreur de type dans __mul__'
+ return t_nappe(self.vale_para,l_fonc,self.para)
+
+ def __repr__(self) :
+ # affichage de la nappe en double colonne
+ texte=[]
+ for i in range(len(self.vale_para)) :
+ texte.append('paramètre : %f' % self.vale_para[i])
+ texte.append(repr(self.l_fonc[i]))
+ return '\n'.join(texte)
+
+ def homo_support(self,other) :
+ # renvoie la nappe self avec un support union de celui de self et de other
+ # le support est la discrétisation vale_para et les discrétisations des fonctions
+ if self==other : return self
+ if self.para!=other.para : raise StandardError, 'combinaison de nappes à caractéristiques interpolation et prolongement différentes'
+ vale_para=self.vale_para.tolist()+other.vale_para.tolist()
+ vale_para=dict([(i,0) for i in vale_para]).keys()
+ vale_para.sort()
+ vale_para=array(vale_para)
+ l_fonc=[]
+ for val in vale_para :
+ if val in self.vale_para : l_fonc.append(self.l_fonc[searchsorted(self.vale_para,val)])
+ elif val in other.vale_para :
+ other_fonc=other.l_fonc[searchsorted(other.vale_para,val)]
+ new_vale_x=other_fonc.vale_x
+ new_para =other_fonc.para
+ new_vale_y=[self(x) for x in new_vale_x]
+ if isinstance(other_fonc,t_fonction) :
+ l_fonc.append(t_fonction(new_vale_x,new_vale_y,new_para))
+ if isinstance(other_fonc,t_fonction_c) :
+ l_fonc.append(t_fonction_c(new_vale_x,new_vale_y,new_para))
+ else : raise StandardError, 'combinaison de nappes : incohérence'
+ return t_nappe(vale_para,l_fonc,self.para)
+
+ def extreme(self) :
+ # renvoie un dictionnaire des valeurs d'ordonnées min et max
+ val_min=min([min(fonc.vale_y) for fonc in self.l_fonc])
+ val_max=max([max(fonc.vale_y) for fonc in self.l_fonc])
+ vm={'min':[],'max':[]}
+ for i in range(len(self.vale_para)) :
+ for j in range(len(self.l_fonc[i].vale_y)) :
+ y = self.l_fonc[i].vale_y[j]
+ if y==val_min : vm['min'].append([y,self.l_fonc[i].vale_x[j],self.vale_para[i]])
+ if y==val_max : vm['max'].append([y,self.l_fonc[i].vale_x[j],self.vale_para[i]])
+ vm['min'].sort()
+ vm['max'].sort()
+ for item in vm['min'] : item.reverse()
+ for item in vm['max'] : item.reverse()
+ return vm
--- /dev/null
+#@ MODIF transpose Utilitai DATE 14/09/2004 AUTEUR MCOURTOI M.COURTOIS
+# -*- coding: iso-8859-1 -*-
+# CONFIGURATION MANAGEMENT OF EDF VERSION
+# ======================================================================
+# COPYRIGHT (C) 1991 - 2004 EDF R&D WWW.CODE-ASTER.ORG
+# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
+# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
+# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
+# (AT YOUR OPTION) ANY LATER VERSION.
+#
+# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
+# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
+# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
+# GENERAL PUBLIC LICENSE FOR MORE DETAILS.
+#
+# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
+# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
+# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
+# ======================================================================
+
+
+
+######################################################################
+#### méthode de transposition de double liste
+#### à résorber quand on aura numarray et possibilité
+#### d opérations sur des arrays de strings
+######################################################################
+def transpose(liste):
+ n=range(len(liste[0]))
+ m=range(len(liste))
+ liste_t=[[] for i in n]
+ for i in n :
+ for j in m :
+ liste_t[i].append(liste[j][i])
+ return liste_t
import ops
# pas d'import aster, Utilitai... pour ne pas gener eficas
-#from Utilitai.t_fonction import t_fonction,t_fonction_c,t_nappe
+from Utilitai.t_fonction import t_fonction,t_fonction_c,t_nappe
try:
import aster
représentation python de la fonction
"""
if arg=='real' :
- from Utilitai.t_fonction import t_fonction,t_fonction_c,t_nappe
return t_fonction(self.Absc(),
self.Ordo(),
self.Parametres())
elif arg=='complex' :
- from Utilitai.t_fonction import t_fonction,t_fonction_c,t_nappe
return t_fonction_c(self.Absc(),
self.Ordo(),
self.Parametres())
Retourne un objet de la classe t_fonction ou t_fonction_c,
représentation python de la fonction complexe
"""
- from Utilitai.t_fonction import t_fonction,t_fonction_c,t_nappe
if arg=='real' :
return t_fonction(self.Absc(),
self.Ordo(),
vale=self.Valeurs()
l_fonc=[]
i=0
- from Utilitai.t_fonction import t_fonction,t_fonction_c,t_nappe
for pf in para[1] :
para_f={'INTERPOL' : pf['INTERPOL_FONC'],
'PROL_DROITE' : pf['PROL_DROITE_FONC'],
),
INFO =SIMP(statut='f',typ='I',defaut= 1,into=( 1 , 2) ),
) ;
-#& MODIF COMMANDE DATE 08/03/2005 AUTEUR LAMARCHE S.LAMARCHE
+#& MODIF COMMANDE DATE 24/05/2005 AUTEUR MABBAS M.ABBAS
# CONFIGURATION MANAGEMENT OF EDF VERSION
# ======================================================================
# COPYRIGHT (C) 1991 - 2001 EDF R&D WWW.CODE-ASTER.ORG
regles=(UN_PARMI('GROUP_MA_ESCL','MAILLE_ESCL'),),
APPARIEMENT =SIMP(statut='f',typ='TXM',defaut="MAIT_ESCL",
into=("NON","NODAL","MAIT_ESCL","MAIT_ESCL_SYME")),
- RECHERCHE =SIMP(statut='f',typ='TXM',defaut="NOEUD_VOISIN",
+ RECHERCHE =SIMP(statut='f',typ='TXM',defaut="NOEUD_BOUCLE",
into=("NOEUD_BOUCLE","NOEUD_VOISIN")),
LISSAGE =SIMP(statut='f',typ='TXM',defaut="NON",
into=("OUI","NON")),
),
INFO =SIMP(statut='f',typ='I',defaut= 1,into=( 1 , 2) ),
) ;
-#& MODIF COMMANDE DATE 28/02/2005 AUTEUR MABBAS M.ABBAS
+#& MODIF COMMANDE DATE 24/05/2005 AUTEUR MABBAS M.ABBAS
# CONFIGURATION MANAGEMENT OF EDF VERSION
# ======================================================================
# COPYRIGHT (C) 1991 - 2001 EDF R&D WWW.CODE-ASTER.ORG
regles=(UN_PARMI('GROUP_MA_ESCL','MAILLE_ESCL'),),
APPARIEMENT =SIMP(statut='f',typ='TXM',defaut="MAIT_ESCL",
into=("NON","NODAL","MAIT_ESCL","MAIT_ESCL_SYME")),
- RECHERCHE =SIMP(statut='f',typ='TXM',defaut="NOEUD_VOISIN",into=("NOEUD_BOUCLE","NOEUD_VOISIN")),
+ RECHERCHE =SIMP(statut='f',typ='TXM',defaut="NOEUD_BOUCLE",into=("NOEUD_BOUCLE","NOEUD_VOISIN")),
LISSAGE =SIMP(statut='f',typ='TXM',defaut="NON",into=("OUI","NON")),
NORMALE =SIMP(statut='f',typ='TXM',defaut="MAIT",into=("MAIT","MAIT_ESCL")),
METHODE =SIMP(statut='f',typ='TXM',defaut="CONTRAINTE",
TITRE =SIMP(statut='f',typ='TXM',max='**'),
INFO =SIMP(statut='f',typ='I',defaut= 1,into=( 1 , 2 ) ),
) ;
-#& MODIF COMMANDE DATE 12/05/2005 AUTEUR DURAND C.DURAND
+#& MODIF COMMANDE DATE 24/05/2005 AUTEUR MCOURTOI M.COURTOIS
# CONFIGURATION MANAGEMENT OF EDF VERSION
# ======================================================================
# COPYRIGHT (C) 1991 - 2001 EDF R&D WWW.CODE-ASTER.ORG
INTERPOL_FONC =SIMP(statut='f',typ='TXM',max=2,into=("NON","LIN","LOG") ),
PROL_DROITE_FONC=SIMP(statut='f',typ='TXM',into=("CONSTANT","LINEAIRE","EXCLU") ),
PROL_GAUCHE_FONC=SIMP(statut='f',typ='TXM',into=("CONSTANT","LINEAIRE","EXCLU") ),
+ INFO =SIMP(statut='f',typ='I',defaut=1,into=(1,2) ),
)
#& MODIF COMMANDE DATE 22/02/2005 AUTEUR DURAND C.DURAND
# CONFIGURATION MANAGEMENT OF EDF VERSION
),
INFO =SIMP(statut='f',typ='I',defaut=1,into=(1,2) ),
) ;
-#& MODIF COMMANDE DATE 22/02/2005 AUTEUR DURAND C.DURAND
+#& MODIF COMMANDE DATE 17/05/2005 AUTEUR CIBHHLV L.VIVAN
# CONFIGURATION MANAGEMENT OF EDF VERSION
# ======================================================================
# COPYRIGHT (C) 1991 - 2001 EDF R&D WWW.CODE-ASTER.ORG
IMPR_GENE=PROC(nom="IMPR_GENE",op= 157,
fr="Calcul du dommage subi par une structure soumise à une sollicitation de type aléatoire",
UIinfo={"groupes":("Impression",)},
+ FORMAT =SIMP(statut='f',typ='TXM',defaut="RESULTAT",into=("RESULTAT",) ),
+ UNITE =SIMP(statut='f',typ='I',defaut=8),
GENE =FACT(statut='o',max='**',
regles=(EXCLUS('TOUT_ORDRE','NUME_ORDRE','INST','FREQ','NUME_MODE',
'LIST_INST','LIST_FREQ','TOUT_MODE','TOUT_INST','LIST_ORDRE'),
EXCLUS('TOUT_PARA','NOM_PARA'),),
# faut-il faire des blocs selon le type de RESU_GENE
RESU_GENE =SIMP(statut='o',typ=(vect_asse_gene_r, tran_gene, mode_gene, harm_gene)),
- FORMAT =SIMP(statut='f',typ='TXM',defaut="RESULTAT",into=("RESULTAT",) ),
- UNITE =SIMP(statut='f',typ='I',defaut=8),
TOUT_ORDRE =SIMP(statut='f',typ='TXM',into=("OUI",) ),
NUME_ORDRE =SIMP(statut='f',typ='I',validators=NoRepeat(),max='**'),
LIST_ORDRE =SIMP(statut='f',typ=listis_sdaster ),
TITRE =SIMP(statut='f',typ='TXM',max='**'),
INFO =SIMP(statut='f',typ='I',defaut=1,into=(1,2) ),
) ;
-#& MODIF COMMANDE DATE 12/05/2005 AUTEUR DURAND C.DURAND
+#& MODIF COMMANDE DATE 24/05/2005 AUTEUR MCOURTOI M.COURTOIS
# CONFIGURATION MANAGEMENT OF EDF VERSION
# ======================================================================
# COPYRIGHT (C) 1991 - 2005 EDF R&D WWW.CODE-ASTER.ORG
CRITERE =SIMP(statut='f',typ='TXM',defaut="RELATIF",into=("RELATIF","ABSOLU") ),
PRECISION =SIMP(statut='f',typ='R',defaut= 1.E-3,val_min=0.E+0 ),
),
+ INFO =SIMP(statut='f',typ='I',defaut=1,into=(1,2) ),
)
#& MODIF COMMANDE DATE 10/06/2004 AUTEUR REZETTE C.REZETTE
# CONFIGURATION MANAGEMENT OF EDF VERSION
-#@ MODIF ops Cata DATE 07/03/2005 AUTEUR DURAND C.DURAND
+#@ MODIF ops Cata DATE 17/05/2005 AUTEUR DURAND C.DURAND
# -*- coding: iso-8859-1 -*-
# CONFIGURATION MANAGEMENT OF EDF VERSION
# ======================================================================
# Modules Eficas
import Accas
from Accas import ASSD
-#from Utilitai.Utmess import UTMESS
+from Utilitai.Utmess import UTMESS
try:
import aster
# On supprime du pickle_context les concepts valant None, ca peut
# etre le cas des concepts non executés, placés après FIN.
pickle_context=get_pickled_context()
- from Utilitai.Utmess import UTMESS
if pickle_context==None :
UTMESS('F','Poursuite',"Erreur a la relecture du fichier pick.1 : aucun objet sauvegardé ne sera récupéré")
return
- from Cata.cata import ASSD
+ from Cata.cata import ASSD,entier
from Noyau.N_CO import CO
for elem in pickle_context.keys():
if type(pickle_context[elem])==types.InstanceType :
if poursu_class!=pickle_class :
UTMESS('F','Poursuite',"Types incompatibles entre glob.1 et pick.1 pour concept de nom "+elem)
return
- elif isinstance(pickle_context[elem],ASSD) and not isinstance(pickle_context[elem],CO) :
+ elif isinstance(pickle_context[elem],ASSD) and pickle_class not in (CO,entier) :
# on n'a pas trouvé le concept dans la base et sa classe est ASSD : ce n'est pas normal
# sauf dans le cas de CO : il n'a alors pas été typé et c'est normal qu'il soit absent de la base
+ # meme situation pour le type 'entier' produit uniquement par DEFI_FICHIER
UTMESS('F','Poursuite',"Concept de nom "+elem+" et de type "+str(pickle_class)+" introuvable dans la base globale")
return
if pickle_context[elem]==None : del pickle_context[elem]
initialdir=os.curdir
# Choix des catalogues
-rep_mat="/local/noyret/Install_Eficas/materiau"
+rep_mat="materiau"
catalogues = (
#('ASTER','v5',os.path.join(rep_cata,'cataSTA5'),'asterv5'),
#('ASTER','v73',os.path.join(rep_cata,'cataSTA73'),'python','defaut'),
('ASTER','v74',os.path.join(rep_cata,'cataSTA74'),'python','defaut'),
('ASTER','v8',os.path.join(rep_cata,'cata_STA8.py'),'python','defaut'),
- ('HOMARD','v1',os.path.join(rep_homard,'homard_cata_V6n.py'),'homard')
)
print "******* Probleme : a la soustraction"
return None
- def __mul__(self,a):
- try :
- return self.valeur*a.valeur
- except :
- print "******* Probleme : a la multiplication"
- return None
-
- def __rmul__(self,a):
- try :
- return self.valeur*a.valeur
- except :
- print "******* Probleme : a la multiplication"
- return None
def __mul__(self,a):
try :
try :
retour = self.valeur * other
except :
- print "******* Probleme : a la multiplication"
+ try :
+ retour = eval(self.valeur) * eval(other)
+ except :
+ try :
+ retour = self.valeur * eval(other)
+ except :
+ print other
+ print "******* Probleme : a la multiplication _mul__"
return retour
def __rmul__ (self,other):
try :
retour = self.valeur * other
except :
- print "******* Probleme : a la multiplication"
+ try :
+ retour = eval(self.valeur) * eval(other)
+ except :
+ print "******* Probleme : a la multiplication __rmul__"
return retour
return retour
+ def cos(self):
+ try :
+ retour=cos(self.valeur)
+ return retour
+ except:
+ print "pb pour cosinus"
+
+ def sin(self):
+ try :
+ retour=sin(self.valeur)
+ return retour
+ except:
+ print "pb pour sinus"
+
+ def tan(self):
+ try :
+ retour=tan(self.valeur)
+ return retour
+ except:
+ print "pb pour tangente"
+
+ def log(self):
+ try :
+ retour=log(self.valeur)
+ return retour
+ except:
+ print "pb pour log"
+
+ def sqrt(self):
+ try :
+ retour=sqrt(self.valeur)
+ return retour
+ except:
+ print "pb pour sqrt"
+
def interprete_valeur(self,val):
"""
Essaie d'interpréter val (chaîne de caractères)comme :
pass
+class COMBI_PARAMETRE :
+ def __init__(self,chainevaleur,valeur):
+ self.chainevaleur=chainevaleur
+ self.valeur=valeur
+
+ def __repr__(self):
+ return self.chainevaleur
+
+ def isvalid(self):
+ if self.valeur and self.chainevaleur:
+ return 1
class ITEM_PARAMETRE :
def __init__(self,param_pere,item=None):
