Aster/Cata/cataSTA76/Macro/ajout_quad_gmsh.py

   1 #@ MODIF ajout_quad_gmsh Macro  DATE 14/09/2004   AUTEUR MCOURTOI M.COURTOIS
   2 # -*- coding: iso-8859-1 -*-
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2002  EDF R&D                  WWW.CODE-ASTER.ORG
   6 # THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
   7 # IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
   8 # THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
   9 # (AT YOUR OPTION) ANY LATER VERSION.
  10 #
  11 # THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
  12 # WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
  13 # MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
  14 # GENERAL PUBLIC LICENSE FOR MORE DETAILS.
  15 #
  16 # YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
  17 # ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
  18 #    1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
  19 # ======================================================================
  20
  21
  22
  23 # script PYTHON de modification d'un fichier de maillage GMSH (*.msh)
  24 # pour transformer les mailles lineiques en mailles quadratiques
  25 # SEG2->SEG3 / TRIA3->TRIA6 / QUAD4->QUAD8,
  26 # TETRA4 -> TETRA10 / PENTA6 -> PENTA15 / PYRAM5 -> PYRAM13
  27
  28 def ajout_quad_gmsh(texte):
  29
  30   import os,sys,copy
  31
  32   try:
  33 # construction du dictionnaire nom du noeud / coordonnees : dict_no
  34 # construction de la liste des noeuds : l_no
  35
  36     typ_mail={}
  37     typ_mail['LI']=['1']
  38     typ_mail['2D']=['2','3']
  39     typ_mail['3D']=['4','5','6','7']
  40     typ_mail['PO']=['15']
  41
  42     texte_eclate=texte.split('\n')
  43     nbno=int(texte_eclate[texte_eclate.index('$NOD')+1])
  44     l_no=texte_eclate[texte_eclate.index('$NOD')+2:texte_eclate.index('$ENDNOD')]
  45     dict_no={}
  46     for i in texte_eclate[texte_eclate.index('$NOD')+2:texte_eclate.index('$ENDNOD')]:
  47       cc=i.split(' ')
  48       dict_no[cc[0]]=[]
  49       for j in cc[1:]:dict_no[cc[0]].append(float(j))
  50
  51     l_el1=texte_eclate[texte_eclate.index('$ELM')+2:texte_eclate.index('$ENDELM')]
  52     nbel =texte_eclate[texte_eclate.index('$ELM')+1]
  53
  54
  55 # construction du dictionnaire : element / liste des aretes de l'element : elems_aretes
  56 # et de son inverse            : aretes / elements contenant cette arete : aretes_elems
  57
  58     aretes_elems={}
  59     elems_aretes={}
  60     l_el=[]
  61     for elem in l_el1 :
  62         connec0=elem.split(' ')[5:]
  63         while '' in connec0 : connec0.remove('')
  64         aa=elem.split(' ')[:4]
  65
  66 # attention : indicateur de type de maille : ajouter 7 pour passer
  67 # de TRIA3 a TRIA6, de SEG2 a SEG3, de QUAD4 a QUAD8 (voir inigms.f)
  68
  69 #   si maille POI1, on ne fait rien
  70
  71         if aa[1] not in typ_mail['PO'] : typel=str(int(aa[1])+7)
  72         else                           : typel=aa[1]
  73         nom_elem=aa[0]+' '+typel+' '+aa[2]+' '+aa[3]
  74         segments=[]
  75
  76 # traitement des mailles lineaires : un seul segment
  77         if   aa[1] in typ_mail['LI']:
  78            segments.append([int(connec0[0]),int(connec0[1])])
  79
  80 # traitement des mailles planes (TRI3 QUA4) : on cree les segments en prenant les noeuds consecutivement, puis on ferme
  81         elif aa[1] in typ_mail['2D']:
  82            for i in range(len(connec0)-1) : segments.append([int(connec0[i]),int(connec0[i+1])])
  83            segments.append([int(connec0[-1]),int(connec0[0])])
  84
  85 # traitement des mailles TETRA4 : pour comprendre, voir "fonctions de forme des elements isoparametriques" R3.01.01
  86         elif aa[1]=='4':
  87            for i in range(0,2) : segments.append([int(connec0[i]),int(connec0[i+1])])
  88            segments.append([int(connec0[2]),int(connec0[0])])
  89            for i in range(0,3) : segments.append([int(connec0[3]),int(connec0[i])])
  90
  91 # traitement des mailles HEXA8 : pour comprendre, voir "fonctions de forme des elements isoparametriques" R3.01.01
  92         elif aa[1]=='5':
  93            for i in range(0,3) : segments.append([int(connec0[i]),int(connec0[i+1])])
  94            segments.append([int(connec0[3]),int(connec0[0])])
  95            for i in range(4,7) : segments.append([int(connec0[i]),int(connec0[i+1])])
  96            segments.append([int(connec0[7]),int(connec0[4])])
  97            for i in range(0,4) : segments.append([int(connec0[i]),int(connec0[i+4])])
  98
  99 # traitement des mailles PENTA6 : pour comprendre, voir "fonctions de forme des elements isoparametriques" R3.01.01
 100         elif aa[1]=='6':
 101            for i in range(0,2) : segments.append([int(connec0[i]),int(connec0[i+1])])
 102            segments.append([int(connec0[2]),int(connec0[0])])
 103            for i in range(3,5) : segments.append([int(connec0[i]),int(connec0[i+1])])
 104            segments.append([int(connec0[5]),int(connec0[3])])
 105            for i in range(0,3) : segments.append([int(connec0[i]),int(connec0[i+3])])
 106
 107 # traitement des mailles PYRA5 : pour comprendre, voir "fonctions de forme des elements isoparametriques" R3.01.01
 108         elif aa[1]=='7':
 109            for i in range(0,3) : segments.append([int(connec0[i]),int(connec0[i+1])])
 110            segments.append([int(connec0[3]),int(connec0[0])])
 111            for i in range(0,4) : segments.append([int(connec0[4]),int(connec0[i])])
 112
 113         elems_aretes[nom_elem]=[]
 114         for seg in segments :
 115             elems_aretes[nom_elem].append(copy.copy(seg))
 116             seg.sort()
 117         for seg in segments :
 118            if aretes_elems.has_key(tuple(seg)):aretes_elems[tuple(seg)].append(nom_elem)
 119            else                               :aretes_elems[tuple(seg)]=[nom_elem]
 120         l_el.append(nom_elem)
 121
 122 # construction de la liste complete des aretes
 123
 124     l_ar=aretes_elems.keys()
 125     l_ar.sort()
 126
 127 # numero du plus grand noeud :
 128     nmax=int(l_no[-1].split(' ')[0])
 129
 130 # nouveau nombre total de noeuds :
 131     ndtot=nbno+len(l_ar)
 132
 133 # insertion des noeuds milieux dans le dictionnaire "aretes" et la liste "l_no"  :
 134     ind=nmax
 135     for i in l_ar:
 136       ind=ind+1
 137       x=(dict_no[str(i[0])][0]+dict_no[str(i[1])][0])/2
 138       y=(dict_no[str(i[0])][1]+dict_no[str(i[1])][1])/2
 139       z=(dict_no[str(i[0])][2]+dict_no[str(i[1])][2])/2
 140       l_no.append(str(ind)+' '+str(x)+' '+str(y)+' '+str(z))
 141       for elem in aretes_elems[i]:
 142         for ar in elems_aretes[elem]:
 143             k=copy.copy(ar)
 144             k.sort()
 145             kk=tuple(k)
 146             if i==kk:
 147               ar.insert(1,ind)
 148
 149 # re-ecriture du fichier avec les noeuds milieux :
 150
 151     resu='$NOD\n'+str(ndtot)+'\n'
 152     for i in l_no:
 153       resu=resu+i+'\n'
 154     resu=resu+'$ENDNOD\n'+'$ELM\n'+nbel+'\n'
 155     for i in l_el:
 156       aa=i.split(' ')[:4]
 157       if aa[1] not in typ_mail['PO']:
 158          typ=str(int(aa[1])-7)
 159       else : typ=aa[1]
 160       if elems_aretes[i]!=[]:
 161         resu=resu+i
 162
 163 # traitement des mailles lineaires : d'abord les noeuds sommets, puis le noeud milieu
 164         if typ in typ_mail['LI']:
 165           resu=resu+' 3 '
 166           for j in elems_aretes[i]:
 167               resu=resu+str(j[0])+' '+str(j[2])+' '+str(j[1])+' '
 168
 169 # traitement des mailles planes (TRI3 QUA4) : d'abord les noeuds sommets, puis, dans le meme ordre, les noeuds milieux
 170         elif typ in typ_mail['2D']:
 171           resu=resu+' '+str(len(elems_aretes[i])*2)+' '
 172           for j in elems_aretes[i]:
 173               resu=resu+str(j[0])+' '
 174           for j in elems_aretes[i]:
 175               resu=resu+str(j[1])+' '
 176
 177 # traitement des mailles TETRA4 : pour comprendre, voir "fonctions de forme des elements isoparametriques" R3.01.01
 178         elif typ=='4':
 179           resu=resu+' 10 '
 180           for j in elems_aretes[i][:4]:
 181               resu=resu+str(j[0])+' '
 182           resu=resu+str(elems_aretes[i][0][1])+' '
 183           resu=resu+str(elems_aretes[i][1][1])+' '
 184           resu=resu+str(elems_aretes[i][2][1])+' '
 185           resu=resu+str(elems_aretes[i][3][1])+' '
 186           resu=resu+str(elems_aretes[i][4][1])+' '
 187           resu=resu+str(elems_aretes[i][5][1])+' '
 188
 189 # traitement des mailles HEXA8 : pour comprendre, voir "fonctions de forme des elements isoparametriques" R3.01.01
 190         elif typ=='5':
 191           resu=resu+' 20 '
 192           for j in elems_aretes[i][:8]:
 193               resu=resu+str(j[0])+' '
 194           resu=resu+str(elems_aretes[i][0][1])+' '
 195           resu=resu+str(elems_aretes[i][1][1])+' '
 196           resu=resu+str(elems_aretes[i][2][1])+' '
 197           resu=resu+str(elems_aretes[i][3][1])+' '
 198
 199           resu=resu+str(elems_aretes[i][8][1])+' '
 200           resu=resu+str(elems_aretes[i][9][1])+' '
 201           resu=resu+str(elems_aretes[i][10][1])+' '
 202           resu=resu+str(elems_aretes[i][11][1])+' '
 203
 204           resu=resu+str(elems_aretes[i][4][1])+' '
 205           resu=resu+str(elems_aretes[i][5][1])+' '
 206           resu=resu+str(elems_aretes[i][6][1])+' '
 207           resu=resu+str(elems_aretes[i][7][1])+' '
 208
 209 # traitement des mailles PENTA6 : pour comprendre, voir "fonctions de forme des elements isoparametriques" R3.01.01
 210         elif typ=='6':
 211           resu=resu+' 15 '
 212           for j in elems_aretes[i][:6]:
 213               resu=resu+str(j[0])+' '
 214           resu=resu+str(elems_aretes[i][0][1])+' '
 215           resu=resu+str(elems_aretes[i][1][1])+' '
 216           resu=resu+str(elems_aretes[i][2][1])+' '
 217
 218           resu=resu+str(elems_aretes[i][6][1])+' '
 219           resu=resu+str(elems_aretes[i][7][1])+' '
 220           resu=resu+str(elems_aretes[i][8][1])+' '
 221
 222           resu=resu+str(elems_aretes[i][3][1])+' '
 223           resu=resu+str(elems_aretes[i][4][1])+' '
 224           resu=resu+str(elems_aretes[i][5][1])+' '
 225
 226 # traitement des mailles PYRA5 : pour comprendre, voir "fonctions de forme des elements isoparametriques" R3.01.01
 227         elif typ=='7':
 228           resu=resu+' 13 '
 229           for j in elems_aretes[i][:4]:
 230               resu=resu+str(j[0])+' '
 231           resu=resu+str(elems_aretes[i][0][1])+' '
 232           resu=resu+str(elems_aretes[i][1][1])+' '
 233           resu=resu+str(elems_aretes[i][2][1])+' '
 234
 235           resu=resu+str(elems_aretes[i][6][1])+' '
 236           resu=resu+str(elems_aretes[i][7][1])+' '
 237           resu=resu+str(elems_aretes[i][8][1])+' '
 238
 239           resu=resu+str(elems_aretes[i][3][1])+' '
 240           resu=resu+str(elems_aretes[i][4][1])+' '
 241           resu=resu+str(elems_aretes[i][5][1])+' '
 242
 243       else:
 244 # traitement des mailles POI1 : on recopie la connectivite telle quelle, sans ajouter de nouveau noeud
 245         resu=resu+l_el1[l_el.index(i)]
 246
 247       resu=resu+'\n'
 248
 249     resu=resu+'$ENDELM\n'
 250     return resu
 251   except :
 252     return 0