[cig-commits] r20756 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES/Gmsh_example_serial UTILS/Gmsh
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Fri Sep 21 05:32:08 PDT 2012
Author: dkomati1
Date: 2012-09-21 05:32:08 -0700 (Fri, 21 Sep 2012)
New Revision: 20756
Added:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/process_the_Gmsh_file_once_and_for_all.sh
Modified:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Surf_abs_SqrCirc
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Surf_free_SqrCirc
seismo/2D/SPECFEM2D/trunk/UTILS/Gmsh/LibGmsh2Specfem_convert_Gmsh_to_Specfem2D_official.py
Log:
updated UTILS/Gmsh/LibGmsh2Specfem_convert_Gmsh_to_Specfem2D_official.py and added EXAMPLES/Gmsh_example_serial/process_the_Gmsh_file_once_and_for_all.sh
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Surf_abs_SqrCirc
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Surf_abs_SqrCirc 2012-09-21 11:49:06 UTC (rev 20755)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Surf_abs_SqrCirc 2012-09-21 12:32:08 UTC (rev 20756)
@@ -1,193 +1,145 @@
-192
-215 2 145 146
-216 2 146 147
-237 2 116 117
-240 2 115 116
-303 2 137 138
-304 2 138 139
-338 2 123 124
-339 2 124 125
-377 2 233 234
-380 2 232 233
-387 2 28 29
-388 2 29 30
-406 2 230 231
-407 2 232 231
-433 2 30 31
-434 2 32 31
-449 2 24 23
-452 2 22 23
-457 2 328 327
-458 2 328 329
-465 2 245 246
-468 2 244 245
-479 2 16 17
-480 2 17 18
-513 2 336 337
-516 2 336 335
-531 2 238 239
-532 2 240 239
-571 2 305 306
-572 2 306 307
-593 2 130 131
-596 2 129 130
-615 2 46 47
-616 2 48 47
-617 2 216 215
-620 2 214 215
-787 2 136 135
-788 2 136 137
-839 2 34 35
-840 2 35 36
-869 2 227 228
-872 2 226 227
-881 2 310 311
-884 2 309 310
-891 2 315 316
-892 2 316 317
-895 2 317 318
-896 2 318 319
-1010 2 321 322
-1011 2 322 323
-1175 2 149 150
-1176 2 150 151
-1189 2 112 113
-1192 2 112 111
-1295 2 331 332
-1296 2 332 333
-1301 2 51 52
-1304 2 50 51
-1315 2 210 211
-1316 2 211 212
-1405 2 41 42
-1408 2 40 41
-1415 2 220 221
-1416 2 221 222
-1453 2 126 127
-1456 2 125 126
-1625 2 21 22
-1628 2 20 21
-1643 2 240 241
-1644 2 241 242
-1646 2 206 207
-1647 2 208 207
-1657 2 54 55
-1658 2 56 55
-1777 2 108 4
-1777 2 59 4
-1778 2 108 109
-1780 2 58 59
-1781 2 154 3
-1781 2 3 203
-1782 2 203 204
-1784 2 153 154
-1786 2 52 53
-1787 2 53 54
-1789 2 208 209
-1790 2 209 210
-1857 2 307 308
-1858 2 308 309
-1886 2 128 127
-1887 2 128 129
-1906 2 333 334
-1907 2 334 335
-1927 2 329 330
-1928 2 330 331
-1965 2 109 110
-1966 2 110 111
-1970 2 152 151
-1971 2 152 153
-2010 2 48 49
-2011 2 49 50
-2013 2 320 321
-2016 2 320 319
-2017 2 212 213
-2018 2 213 214
-2109 2 340 341
-2112 2 339 340
-2115 2 301 302
-2116 2 302 303
-2154 2 113 114
-2155 2 114 115
-2301 2 38 39
-2302 2 40 39
-2304 2 222 223
-2305 2 224 223
-2346 2 133 134
-2347 2 134 135
-2379 2 26 27
-2380 2 27 28
-2382 2 234 235
-2383 2 235 236
-2430 2 218 219
-2431 2 219 220
-2433 2 42 43
-2434 2 43 44
-2448 2 120 119
-2449 2 120 121
-2451 2 141 142
-2452 2 142 143
-2490 2 117 118
-2491 2 118 119
-2493 2 144 143
-2494 2 144 145
-2517 2 325 326
-2518 2 326 327
-2520 2 337 338
-2521 2 338 339
-2523 2 304 303
-2524 2 304 305
-2562 2 236 237
-2563 2 237 238
-2565 2 24 25
-2566 2 25 26
-2568 2 131 132
-2569 2 132 133
-2586 2 216 217
-2587 2 217 218
-2589 2 44 45
-2590 2 45 46
-2592 2 312 311
-2593 2 312 313
-2610 2 323 324
-2611 2 324 325
-2622 2 139 140
-2623 2 140 141
-2625 2 121 122
-2626 2 122 123
-2628 2 14 15
-2629 2 16 15
-2631 2 246 247
-2632 2 248 247
-2646 2 313 314
-2647 2 314 315
-2649 2 242 243
-2650 2 243 244
-2652 2 18 19
-2653 2 19 20
-2661 2 32 33
-2662 2 33 34
-2664 2 228 229
-2665 2 229 230
-2673 2 147 148
-2674 2 148 149
-2676 2 36 37
-2677 2 37 38
-2679 2 224 225
-2680 2 225 226
-2700 2 344 1
-2702 2 344 343
-2703 2 249 2
-2705 2 248 249
-2706 2 1 13
-2707 2 13 14
-2709 2 2 298
-2710 2 298 299
-2712 2 56 57
-2713 2 57 58
-2715 2 204 205
-2716 2 205 206
-2727 2 341 342
-2728 2 342 343
-2730 2 299 300
-2731 2 300 301
+144
+215 2 145 146 4
+216 2 146 147 4
+237 2 116 117 4
+240 2 115 116 4
+303 2 137 138 4
+304 2 138 139 4
+338 2 123 124 4
+339 2 124 125 4
+377 2 233 234 1
+380 2 232 233 1
+406 2 230 231 1
+407 2 232 231 1
+457 2 328 327 2
+458 2 328 329 2
+465 2 245 246 1
+468 2 244 245 1
+513 2 336 337 2
+516 2 336 335 2
+531 2 238 239 1
+532 2 240 239 1
+571 2 305 306 2
+572 2 306 307 2
+593 2 130 131 4
+596 2 129 130 4
+617 2 216 215 1
+620 2 214 215 1
+787 2 136 135 4
+788 2 136 137 4
+869 2 227 228 1
+872 2 226 227 1
+881 2 310 311 2
+884 2 309 310 2
+891 2 315 316 2
+892 2 316 317 2
+895 2 317 318 2
+896 2 318 319 2
+1010 2 321 322 2
+1011 2 322 323 2
+1175 2 149 150 4
+1176 2 150 151 4
+1189 2 112 113 4
+1192 2 112 111 4
+1295 2 331 332 2
+1296 2 332 333 2
+1315 2 210 211 1
+1316 2 211 212 1
+1415 2 220 221 1
+1416 2 221 222 1
+1453 2 126 127 4
+1456 2 125 126 4
+1643 2 240 241 1
+1644 2 241 242 1
+1646 2 206 207 1
+1647 2 208 207 1
+1777 2 108 4 4
+1778 2 108 109 4
+1781 2 3 203 1
+1781 2 154 3 4
+1782 2 203 204 1
+1784 2 153 154 4
+1789 2 208 209 1
+1790 2 209 210 1
+1857 2 307 308 2
+1858 2 308 309 2
+1886 2 128 127 4
+1887 2 128 129 4
+1906 2 333 334 2
+1907 2 334 335 2
+1927 2 329 330 2
+1928 2 330 331 2
+1965 2 109 110 4
+1966 2 110 111 4
+1970 2 152 151 4
+1971 2 152 153 4
+2013 2 320 321 2
+2016 2 320 319 2
+2017 2 212 213 1
+2018 2 213 214 1
+2109 2 340 341 2
+2112 2 339 340 2
+2115 2 301 302 2
+2116 2 302 303 2
+2154 2 113 114 4
+2155 2 114 115 4
+2304 2 222 223 1
+2305 2 224 223 1
+2346 2 133 134 4
+2347 2 134 135 4
+2382 2 234 235 1
+2383 2 235 236 1
+2430 2 218 219 1
+2431 2 219 220 1
+2448 2 120 119 4
+2449 2 120 121 4
+2451 2 141 142 4
+2452 2 142 143 4
+2490 2 117 118 4
+2491 2 118 119 4
+2493 2 144 143 4
+2494 2 144 145 4
+2517 2 325 326 2
+2518 2 326 327 2
+2520 2 337 338 2
+2521 2 338 339 2
+2523 2 304 303 2
+2524 2 304 305 2
+2562 2 236 237 1
+2563 2 237 238 1
+2568 2 131 132 4
+2569 2 132 133 4
+2586 2 216 217 1
+2587 2 217 218 1
+2592 2 312 311 2
+2593 2 312 313 2
+2610 2 323 324 2
+2611 2 324 325 2
+2622 2 139 140 4
+2623 2 140 141 4
+2625 2 121 122 4
+2626 2 122 123 4
+2631 2 246 247 1
+2632 2 248 247 1
+2646 2 313 314 2
+2647 2 314 315 2
+2649 2 242 243 1
+2650 2 243 244 1
+2664 2 228 229 1
+2665 2 229 230 1
+2673 2 147 148 4
+2674 2 148 149 4
+2679 2 224 225 1
+2680 2 225 226 1
+2700 2 344 1 2
+2702 2 344 343 2
+2703 2 249 2 1
+2705 2 248 249 1
+2709 2 2 298 2
+2710 2 298 299 2
+2715 2 204 205 1
+2716 2 205 206 1
+2727 2 341 342 2
+2728 2 342 343 2
+2730 2 299 300 2
+2731 2 300 301 2
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Surf_free_SqrCirc
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Surf_free_SqrCirc 2012-09-21 11:49:06 UTC (rev 20755)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/Surf_free_SqrCirc 2012-09-21 12:32:08 UTC (rev 20756)
@@ -1 +1,49 @@
-0
+48
+387 2 28 29
+388 2 29 30
+433 2 30 31
+434 2 32 31
+449 2 24 23
+452 2 22 23
+479 2 16 17
+480 2 17 18
+615 2 46 47
+616 2 48 47
+839 2 34 35
+840 2 35 36
+1301 2 51 52
+1304 2 50 51
+1405 2 41 42
+1408 2 40 41
+1625 2 21 22
+1628 2 20 21
+1657 2 54 55
+1658 2 56 55
+1777 2 59 4
+1780 2 58 59
+1786 2 52 53
+1787 2 53 54
+2010 2 48 49
+2011 2 49 50
+2301 2 38 39
+2302 2 40 39
+2379 2 26 27
+2380 2 27 28
+2433 2 42 43
+2434 2 43 44
+2565 2 24 25
+2566 2 25 26
+2589 2 44 45
+2590 2 45 46
+2628 2 14 15
+2629 2 16 15
+2652 2 18 19
+2653 2 19 20
+2661 2 32 33
+2662 2 33 34
+2676 2 36 37
+2677 2 37 38
+2706 2 1 13
+2707 2 13 14
+2712 2 56 57
+2713 2 57 58
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/process_the_Gmsh_file_once_and_for_all.sh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/process_the_Gmsh_file_once_and_for_all.sh (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/process_the_Gmsh_file_once_and_for_all.sh 2012-09-21 12:32:08 UTC (rev 20756)
@@ -0,0 +1,7 @@
+#!/bin/bash
+#
+# create the absorbing and free surface files from the Gmsh file
+#
+
+python ../../UTILS/Gmsh/LibGmsh2Specfem_convert_Gmsh_to_Specfem2D_official.py SqrCirc.msh -t F -l A -b A -r A
+
Property changes on: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Gmsh_example_serial/process_the_Gmsh_file_once_and_for_all.sh
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/2D/SPECFEM2D/trunk/UTILS/Gmsh/LibGmsh2Specfem_convert_Gmsh_to_Specfem2D_official.py
===================================================================
--- seismo/2D/SPECFEM2D/trunk/UTILS/Gmsh/LibGmsh2Specfem_convert_Gmsh_to_Specfem2D_official.py 2012-09-21 11:49:06 UTC (rev 20755)
+++ seismo/2D/SPECFEM2D/trunk/UTILS/Gmsh/LibGmsh2Specfem_convert_Gmsh_to_Specfem2D_official.py 2012-09-21 12:32:08 UTC (rev 20756)
@@ -1,15 +1,18 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
-# Python code to link gmsh with specfem
+# Python code to link gmsh with specfem
+# New version with identification of the zone
+# and handling of PML
#
#@author: Cristini Paul,
# Laboratoire de Mecanique et d'Acoustique, CNRS, Marseille, France
#
-# February 2012
+# September 2012, 6rd
#
-import sys, string, time
-from os.path import splitext, isfile
+import sys, string
+from os.path import splitext
+
try:
from numpy import *
except ImportError:
@@ -21,37 +24,41 @@
# Ct le nombre de lignes a lire
# Var est le nom de la variable contenant les informations a ecrire
# Fv est le format d'ecriture '%f' pour les noeuds, '%i' pour les autres fichiers
- savetxt(Ng,(Ct,), fmt='%i')
+ savetxt(Ng,(Ct,), fmt='%i') # on met le nombre de lignes en entete
fd = open(Ng,'a')
savetxt(fd, Var, fmt=Fv)
fd.close()
return
#
def OuvreGmsh(Dir,Nom,Bords):
- # Lecture de fichiers .msh genere avec Gmsh
+ # Reading file .msh created with Gmsh
if splitext(Nom)[-1]=='.msh':
fic=Nom
+ Nom=splitext(Nom)[0]
elif splitext(Nom)[-1]=='':
fic=Nom+'.msh'
else:
print 'File extension is not correct'
print 'script aborted'
sys.exit()
+ print '#'*20+' STARTING PROCESSING '+'#'*20
+ print 'File : ', fic
+ print '-'*60
#
# Open the file and get the lines
#
f = file(Dir+fic,'r')
lignes= f.readlines()
f.close()
- # Recherche des positions
- #MotsCles=['','']
+ # Looking for positions
+ #
for ii in range(len(lignes)):
if lignes[ii]=='$Nodes\n': PosNodes=ii
if lignes[ii]=='$PhysicalNames\n': PosPhys=ii
if lignes[ii]=='$Elements\n':
PosElem=ii
break
- # Type d'elements 4 noeuds ou 9 noeuds
+ # Elements type 4 nodes or 9 nodes
TypElem1D = int(string.split(lignes[PosElem+2])[1])
if TypElem1D==1:
Ngnod, LinElem, SurfElem = 4, 1, 3
@@ -60,307 +67,309 @@
Ngnod, LinElem, SurfElem = 9, 8, 10
len1D, len2D = 3, 9
else:
- print 'Element type is not 4 nor 9 nodes'
+ print 'The number of nodes of elemnts is not 4 nor 9'
#-------------------------------------------------------------------------
- # Conditions aux bords du domaine
- # Possible choices: Abso, Free or Perio
- Bord_abso, Bord_free = [], [] # Initialisation
+ # Conditions on sides of the domain
+ # Possible choices: Abso, Free
+ Bord_abso, Bord_free = [], [] # Initialization
+ PML = False
print Bords
#-------------------------------------------------------------------------
# PHYSICAL NAMES
NbPhysNames = int(string.split(lignes[PosPhys+1])[0])
- print 'PhysNames', NbPhysNames
+ # Detection of the presence of PML
+ for Ip in range(NbPhysNames):
+ Nam = string.split(lignes[PosPhys+2+Ip])[2][1:-1]
+ if Nam.endswith('PML'):
+ PML = True
+ break
+ print 'Number of physical Names : ', NbPhysNames
+ print 'PML : ', PML
# Structure de la variable : 1 entier et une chaine de caractere de taille 16
dt = dtype([('dimension',int), ('zone', int), ('name', str, 16)])
PhysCar=zeros((NbPhysNames,), dtype=dt)
+ Med=zeros(NbPhysNames+1) # Les milieux: domaines dont le nom commence par M
+ #--------------------------------------------------------------------------
+ # Initialisation des zones (celles qui ne seront pas affectées auront
+ # une valeur nulle)
+ PML_right, PML_right_bottom, PML_right_top = 0, 0, 0
+ PML_left, PML_left_bottom, PML_left_top = 0, 0, 0
+ PML_top, PML_bottom = 0, 0
+
for Ip in range(NbPhysNames):
Dim = int(string.split(lignes[PosPhys+2+Ip])[0])
Zon = int(string.split(lignes[PosPhys+2+Ip])[1])
Nam = string.split(lignes[PosPhys+2+Ip])[2][1:-1]
PhysCar[Ip] = (Dim, Zon, Nam)
+ print PhysCar[Ip]
+ if Nam.startswith('M'): Med[Zon]=int(Nam[1:]) # Milieux
if Bords.has_key(Nam):
if Bords[Nam] == 'Abso': Bord_abso.append(Zon)
if Bords[Nam] == 'Free': Bord_free.append(Zon)
- if Nam == 'Right': Bord_right=Zon
- if Nam == 'Left': Bord_left=Zon
- if Nam == 'Top': Bord_top=Zon
+ if PML:
+ #------------------------------------------------------------------
+ # Les bords du domaine PML
+ if Nam == 'Right_PML' : Bord_right=Zon
+ if Nam == 'Left_PML' : Bord_left=Zon
+ if Nam == 'Top_PML' : Bord_top=Zon
+ if Nam == 'Bottom_PML': Bord_bottom=Zon
+ else:
+ #------------------------------------------------------------------
+ # Les bords du domaine initial
+ if Nam == 'Right' : Bord_right=Zon
+ if Nam == 'Left' : Bord_left=Zon
+ if Nam == 'Top' : Bord_top=Zon
if Nam == 'Bottom': Bord_bottom=Zon
- #--------------------------------------
- print 'Physical Names', PhysCar
- print 'Absorbing boundaries', Bord_abso
- print 'Free boundaries', Bord_free
- print 'Right boundaries', Bord_right
- print 'Left boundaries', Bord_left
- print 'Top boundaries', Bord_top
- print 'Bottom boundaries', Bord_bottom
+ if PML:
+ #------------------------------------------------------------------
+ # Les differents domaines PML
+ if Nam == 'R' : PML_right, Med[Zon] = Zon, 1
+ if Nam == 'RB' : PML_right_bottom, Med[Zon] = Zon, 1 # Corner
+ if Nam == 'L' : PML_left, Med[Zon] = Zon, 1
+ if Nam == 'LB' : PML_left_bottom, Med[Zon] = Zon, 1 # Corner
+ if Nam == 'T' : PML_top, Med[Zon] = Zon, 1
+ if Nam == 'RT' : PML_right_top, Med[Zon] = Zon, 1 # Corner
+ if Nam == 'LT' : PML_left_top, Med[Zon] = Zon, 1 # Corner
+ if Nam == 'B' : PML_bottom, Med[Zon] = Zon, 1
+ #--------------------------------------------------------------------------
+ print '-'*60
+ print 'Absorbing boundaries : ', Bord_abso
+ print 'Free boundaries : ', Bord_free
+ print '-'*60
+ print 'Right boundaries : ', Bord_right
+ print 'Left boundaries : ', Bord_left
+ print 'Top boundaries : ', Bord_top
+ print 'Bottom boundaries : ', Bord_bottom
+ print '-'*60
+ print 'Med : ',Med
#---------------------------------------------------------------------------
# Infos sur le fichier Gmsh
Ver=float(string.split(lignes[1])[0])
File_Type=int(string.split(lignes[1])[1])
Data_Size=int(string.split(lignes[1])[2])
- # Lecture de noeuds
+ # Nodes
NbNodes=int(string.split(lignes[PosNodes+1])[0])
print 'Number of nodes: ',NbNodes
Nodes=zeros((NbNodes,2),dtype=float)
for Ninc in range(NbNodes):
- Nodes[Ninc]= [float(val) for val in (string.split(lignes[PosNodes+2+Ninc])[1:3])]
+ Nodes[Ninc]= [float(val) for val in \
+ (string.split(lignes[PosNodes+2+Ninc])[1:3])]
#
- # Sauvegarde au format SPECFEM
+ # Save to SPECFEM format
SauvFicSpecfem('Nodes_'+Nom, NbNodes, Nodes, '%f')
- # Lecture des elements
+ # Elements
DecElem=12+NbNodes
NbElements=int(string.split(lignes[PosElem+1])[0])
- print 'Number of elements: ', NbElements
- # depend de l'ordre
- Elements = empty((NbElements,len2D),dtype=int)
- Milieu = empty((NbElements,1),dtype=int)
- Elements2DBord = empty((NbElements),dtype=int)
- Elements1D = empty((NbElements,len1D),dtype=int)
- Elements1DBord = empty((NbElements,len1D),dtype=int)
- #---------------------------------------------------------------------------
+ print 'Number of elements: ', NbElements
+ #--------------------------------------------------------------------------
+ # Initialization
+ #
+ Elements = empty((NbElements,len2D),dtype=int)
+ Milieu = empty((NbElements,1),dtype=int)
+ Elements1D = empty((NbElements,len1D),dtype=int)
+ #--------------------------------------------------------------------------
Elements1DBordAbso = empty((NbElements,len1D),dtype=int)
Elements1DBordFree = empty((NbElements,len1D),dtype=int)
- Elements2DBordAbso = zeros((NbElements,4),dtype=int)
+ Elements2DBordAbso = zeros((NbElements,5),dtype=int) # only 2 nodes + zone
Elements2DBordFree = zeros((NbElements,4),dtype=int)
- #---------------------------------------------------------------------------
- Elements1DBordTop = empty((NbElements,len1D),dtype=int)
+ #--------------------------------------------------------------------------
+ Elements1DBordTop = empty((NbElements,len1D),dtype=int)
Elements1DBordBottom = empty((NbElements,len1D),dtype=int)
- Elements2DBordTop = zeros((NbElements,4),dtype=int)
- Elements2DBordBottom = zeros((NbElements,4),dtype=int)
- #---------------------------------------------------------------------------
+ #--------------------------------------------------------------------------
Elements1DBordRight = empty((NbElements,len1D),dtype=int)
- Elements1DBordLeft = empty((NbElements,len1D),dtype=int)
- Elements2DBordRight = zeros((NbElements,4),dtype=int)
- Elements2DBordLeft = zeros((NbElements,4),dtype=int)
- #---------------------------------------------------------------------------
- DecElem+=1
- Ninc1D, Ninc2D, Ninc1DBord, Ninc1DBordAbso, Ninc1DBordFree = 0, 0, 0, 0, 0
- Ninc1DBordLeft, Ninc1DBordRight, Ninc1DBordTop, Ninc1DBordBottom = 0, 0, 0, 0
+ Elements1DBordLeft = empty((NbElements,len1D),dtype=int)
+ #--------------------------------------------------------------------------
+ if PML:
+ ElementsPML = empty((NbElements,2),dtype=int) # éléments + codage
+ DecElem+=1 #
+ #--------------------------------------------------------------------------
+ # Initialization of counters
+ N1D, N2D, N2DPML, N1DBord, N1DBordAbso, N1DBordFree = 0, 0, 0, 0, 0 ,0
+ N1DBordLeft, N1DBordRight, N1DBordTop, N1DBordBottom = 0, 0, 0, 0
+ #--------------------------------------------------------------------------
+ # Loop over all elements
+ #
for Ninc in range(NbElements):
- Pos = PosElem+Ninc+2
+ Pos = PosElem+Ninc+2
+ Ispec = int(string.split(lignes[Pos])[0])
TypElem = int(string.split(lignes[Pos])[1])
ZonP = int(string.split(lignes[Pos])[3])
- Milieu[Ninc2D]= 1
- if TypElem==LinElem:
- Elements1D[Ninc1D] = [int(val) for val in (string.split(lignes[Pos])[5:])]
- # Bord droit
+ #
+ if TypElem==LinElem: # Elements 1D (lines)
+ Elements1D[N1D] = [int(val) for val in \
+ (string.split(lignes[Pos])[5:])]
+ # Bottom
+ if ZonP==Bord_bottom:
+ Elements1DBordBottom[N1DBordBottom] = Elements1D[N1D]
+ N1DBordBottom+=1
+ # Top
+ if ZonP==Bord_top:
+ Elements1DBordTop[N1DBordTop] = Elements1D[N1D]
+ N1DBordTop+=1
+ # Left
+ if ZonP==Bord_left:
+ Elements1DBordLeft[N1DBordLeft] = Elements1D[N1D]
+ N1DBordLeft+=1
+ # Right
if ZonP==Bord_right:
- Elements1DBordRight[Ninc1DBordRight] = Elements1D[Ninc1D]
- Ninc1DBordRight+=1
- Elements1DBord[Ninc1DBord] = Elements1D[Ninc1D]
- Ninc1DBord+=1
+ Elements1DBordRight[N1DBordRight] = Elements1D[N1D]
+ N1DBordRight+=1
+ N1D+=1
+
+ #------------------------------------(---------------------------------
+ if TypElem==SurfElem:
+ Elements[N2D]= [int(val) for val in \
+ (string.split(lignes[Pos])[5:])]
+ Milieu[N2D]= Med[ZonP]
+ if PML:
+ # PML Bottom
+ if ZonP==PML_bottom:
+ ElementsPML[N2DPML,:] = [N2D+1, 8]
+ N2DPML+=1
+ # PML Right
+ if ZonP==PML_right:
+ ElementsPML[N2DPML,:] = [N2D+1, 2]
+ N2DPML+=1
+ # PML Top
+ if ZonP==PML_top:
+ ElementsPML[N2DPML,:] = [N2D+1, 4]
+ N2DPML+=1
+ # PML Left
+ if ZonP==PML_left:
+ ElementsPML[N2DPML,:] = [N2D+1, 1]
+ N2DPML+=1
+ #--------------------------------------------------------------
+ # PML Corner Right Bottom
+ if ZonP==PML_right_bottom:
+ ElementsPML[N2DPML,:] = [N2D+1, 10]
+ N2DPML+=1
+ # PML Corner Right Top
+ if ZonP==PML_right_top:
+ ElementsPML[N2DPML,:] = [N2D+1, 6]
+ N2DPML+=1
+ # PML Corner Left Top
+ if ZonP==PML_left_top:
+ ElementsPML[N2DPML,:] = [N2D+1, 5]
+ N2DPML+=1
+ # PML Corner Left Bottom
+ if ZonP==PML_left_bottom:
+ ElementsPML[N2DPML,:] = [N2D+1, 9]
+ N2DPML+=1
+ N2D+=1
+ #--------------------------------------------------------------------------
+ print '-'*60
+ print "Elements 1D, 2D : ", N1D, N2D
+ if PML: print "Elements PML : ", N2DPML
+ #--------------------------------------------------------------------------
+ Elements = Elements[:N2D,:]
+ Milieu = Milieu[:N2D,:]
+ if PML: ElementsPML = ElementsPML[:N2DPML,:]
+ # removal of unnecessary zeros created at initialization
+ Elements1D = Elements1D[:N1D,:]
+ #
+ Elements1DBordAbso = Elements1DBordAbso[:N1DBordAbso,:]
+ Elements1DBordFree = Elements1DBordFree[:N1DBordFree,:]
+ #
+ Elements1DBordRight = Elements1DBordRight[:N1DBordRight,:]
+ Elements1DBordLeft = Elements1DBordLeft[:N1DBordLeft,:]
+ Elements1DBordTop = Elements1DBordTop[:N1DBordTop,:]
+ Elements1DBordBottom = Elements1DBordBottom[:N1DBordBottom,:]
+
+ # Creation of a vector containing only the nodes of one side
+ # We take the first two nodes (end nodes) not intermediate (high order)
+ #--------------------------------------------------------------------------
+ NodesBordRight = set(ravel(Elements1DBordRight[:,:2]))
+ NodesBordLeft = set(ravel(Elements1DBordLeft[:,:2]))
+ NodesBordTop = set(ravel(Elements1DBordTop[:,:2]))
+ NodesBordBottom = set(ravel(Elements1DBordBottom[:,:2]))
+ #--------------------------------------------------------------------------
+ ctf, cta = 0, 0
+ #
+ print '#'*60
+ print 'Identification of elements in contact with sides'
+ print '# Corners '
+ #--------------------------------------------------------------------------
+ # Detection of elements in contact with sides
+
+ for Ct2D in xrange(N2D):
+ #
+ jj=set(Elements[Ct2D,:]) # take the nodes of the element
+ ct_corner = 0
+
+ # Bottom
+ if not set.isdisjoint(jj, NodesBordBottom):
+ rr = set.intersection(jj, NodesBordBottom)
+ lrr = len(rr)
+ if lrr > 1:
+ ela = concatenate(([Ct2D+1], [lrr], list(rr), [1]))
+ elf = concatenate(([Ct2D+1], [lrr], list(rr)))
+ if Bords['Bottom']=='Abso':
+ Elements2DBordAbso[cta,:] = ela
+ cta+=1
+ if Bords['Bottom']=='Free':
+ Elements2DBordFree[ctf,:] = elf
+ ctf+=1
+ ct_corner +=1
+
+ # Right
+ if not set.isdisjoint(jj, NodesBordRight):
+ rr = set.intersection(jj, NodesBordRight)
+ lrr = len(rr)
+ if lrr > 1:
+ ela = concatenate(([Ct2D+1], [lrr], list(rr), [2]))
+ elf = concatenate(([Ct2D+1], [lrr], list(rr)))
if Bords['Right']=='Abso':
- Elements1DBordAbso[Ninc1DBordAbso] = Elements1D[Ninc1D]
- Ninc1DBordAbso+=1
- else:
- Elements1DBordFree[Ninc1DBordFree] = Elements1D[Ninc1D]
- Ninc1DBordFree+=1
- # Bord gauche
- if ZonP==Bord_left:
- Elements1DBordLeft[Ninc1DBordLeft] = Elements1D[Ninc1D]
- Ninc1DBordLeft+=1
- Elements1DBord[Ninc1DBord] = Elements1D[Ninc1D]
- Ninc1DBord+=1
+ Elements2DBordAbso[cta,:] = ela
+ cta+=1
+ if Bords['Right']=='Free':
+ Elements2DBordFree[ctf,:] = elf
+ ctf+=1
+ ct_corner +=1
+
+ # Top
+ if not set.isdisjoint(jj, NodesBordTop):
+ rr = set.intersection(jj, NodesBordTop)
+ lrr = len(rr)
+ if lrr > 1:
+ ela = concatenate(([Ct2D+1], [lrr], list(rr), [3]))
+ elf = concatenate(([Ct2D+1], [lrr], list(rr)))
+ if Bords['Top']=='Abso':
+ Elements2DBordAbso[cta,:] = ela
+ cta+=1
+ if Bords['Top']=='Free':
+ Elements2DBordFree[ctf,:] = elf
+ ctf+=1
+ ct_corner +=1
+
+ # Left
+ if not set.isdisjoint(jj, NodesBordLeft):
+ rr = set.intersection(jj, NodesBordLeft)
+ lrr = len(rr)
+ if lrr > 1:
+ ela = concatenate(([Ct2D+1], [lrr], list(rr), [4]))
+ elf = concatenate(([Ct2D+1], [lrr], list(rr)))
if Bords['Left']=='Abso':
- Elements1DBordAbso[Ninc1DBordAbso] = Elements1D[Ninc1D]
- Ninc1DBordAbso+=1
- else:
- Elements1DBordFree[Ninc1DBordFree] = Elements1D[Ninc1D]
- Ninc1DBordFree+=1
- # Bord haut
- if ZonP==Bord_top:
- Elements1DBordTop[Ninc1DBordTop] = Elements1D[Ninc1D]
- Ninc1DBordTop+=1
- Elements1DBord[Ninc1DBord] = Elements1D[Ninc1D]
- Ninc1DBord+=1
- if Bords['Top']=='Abso':
- Elements1DBordAbso[Ninc1DBordAbso] = Elements1D[Ninc1D]
- Ninc1DBordAbso+=1
- else:
- Elements1DBordFree[Ninc1DBordFree] = Elements1D[Ninc1D]
- Ninc1DBordFree+=1
- # Bord bas
- if ZonP==Bord_bottom:
- Elements1DBordBottom[Ninc1DBordBottom] = Elements1D[Ninc1D]
- Ninc1DBordBottom+=1
- Elements1DBord[Ninc1DBord] = Elements1D[Ninc1D]
- Ninc1DBord+=1
- if Bords['Bottom']=='Abso':
- Elements1DBordAbso[Ninc1DBordAbso] = Elements1D[Ninc1D]
- Ninc1DBordAbso+=1
- else:
- Elements1DBordFree[Ninc1DBordFree] = Elements1D[Ninc1D]
- Ninc1DBordFree+=1
- Ninc1D+=1
- if TypElem==SurfElem:
- Elements[Ninc2D]= [int(val) for val in (string.split(lignes[Pos])[5:])]
- Milieu[Ninc2D]= ZonP-4
- Ninc2D+=1
- Elements = Elements[:Ninc2D,:]
- Milieu = Milieu[:Ninc2D,:]
+ Elements2DBordAbso[cta,:] = ela
+ cta+=1
+ if Bords['Left']=='Free':
+ Elements2DBordFree[ctf,:] = elf
+ ctf+=1
+ ct_corner +=1
+
+ if ct_corner > 1: print 'Corner', ct_corner, Ct2D+1
+
#
- Elements1D=Elements1D[:Ninc1D,:]
- Elements1DBord=Elements1DBord[:Ninc1DBord,:]
- Elements1DBordAbso=Elements1DBordAbso[:Ninc1DBordAbso,:]
- Elements1DBordFree=Elements1DBordFree[:Ninc1DBordFree,:]
- Elements1DBordRight=Elements1DBordRight[:Ninc1DBordRight,:]
- Elements1DBordLeft=Elements1DBordLeft[:Ninc1DBordLeft,:]
- Elements1DBordTop=Elements1DBordTop[:Ninc1DBordTop,:]
- Elements1DBordBottom=Elements1DBordBottom[:Ninc1DBordBottom,:]
- # Modification la matrice contenant les elements 1D pour travailler sur les noeuds
- Elements1DBordFlat=ravel(Elements1DBord)
- # Noeuds appartenant aux bords du domaine
- NodesBordC=set(Elements1DBordFlat) # permet d'enlever les elements dupliques
- NodesBord2n=set(Elements1DBordFlat[::3]) # Noeuds aux bords sans intermediaire
- #-------------------------------------------------------
- NodesBordRight = set(ravel(Elements1DBordRight))
- NodesBordLeft = set(ravel(Elements1DBordLeft))
- NodesBordTop = set(ravel(Elements1DBordTop))
- NodesBordBottom = set(ravel(Elements1DBordBottom))
- #-------------------------------------------------------
- Elt=ravel(Elements1DBordRight)
- NodesBordRight2n = set(hstack((Elt,Elt[-2]))[::3])
- Elt=ravel(Elements1DBordLeft)
- NodesBordLeft2n = set(hstack((Elt,Elt[-2]))[::3])
- Elt=ravel(Elements1DBordTop)
- NodesBordTop2n = set(hstack((Elt,Elt[-2]))[::3])
- Elt=ravel(Elements1DBordBottom)
- NodesBordBottom2n = set(hstack((Elt,Elt[-2]))[::3])
- #
- ctBord=0
- # Detection des elements qui sont aux bords droit gauche haut bas
- t0=time.clock()
- ctf, cta, ctr, ctl, ctt, ctb = 0, 0, 0, 0, 0, 0
- for Ct2D in xrange(Ninc2D):
- # Attention aux elements qui ont un point correspondant a un coin
- jj=set(Elements[Ct2D])
- if not set.isdisjoint(jj, NodesBordC):
- ctBord+=1
- # Bord Right
- if not set.isdisjoint(jj, NodesBordRight2n):
- rr = set.intersection(jj, NodesBordRight2n)
- if len(rr)==2:
- el = concatenate(([Ct2D+1], [len(rr)], list(rr)))
- Elements2DBordRight[ctr,:] = el
- ctr+=1
- if Bords['Right']=='Abso':
- Elements2DBordAbso[cta,:] = el
- cta+=1
- if Bords['Right']=='Free':
- Elements2DBordFree[ctf,:] = el
- ctf+=1
- # Cas du bord double
- if len(rr)==3:
- ctdb=0
- # Recherche des deux elements 1D
- for Nn in xrange(Ninc1DBordRight):
- kk = set(Elements1DBordRight[Nn])
- if not set.isdisjoint(rr, kk):
- if kk.issubset(jj):
- exec 'kk'+str(ctdb)+'=kk'
- ctdb+=1
- Elc = set.intersection(kk0, kk1)
- rrm=rr.difference(Elc)
- rr1=list(rrm)
- El1=concatenate(([rr1[0]],list(Elc)))
- El2=concatenate((list(Elc),[rr1[1]]))
- el1 = concatenate(([Ct2D+1], [2], El1))
- el2 = concatenate(([Ct2D+1], [2], El2))
- Elements2DBordRight[ctr,:] = el1
- ctr+=1
- Elements2DBordRight[ctr,:] = el2
- ctr+=1
- if Bords['Right']=='Abso':
- Elements2DBordAbso[cta,:] = el1
- cta+=1
- Elements2DBordAbso[cta,:] = el2
- cta+=1
- if Bords['Right']=='Free':
- Elements2DBordFree[ctf,:] = el1
- ctf+=1
- Elements2DBordFree[ctf,:] = el2
- ctf+=1
- # Bord Left
- if not set.isdisjoint(jj, NodesBordLeft2n):
- rr = set.intersection(jj, NodesBordLeft2n)
- if len(rr)==2:
- el = concatenate(([Ct2D+1], [len(rr)], list(rr)))
- Elements2DBordLeft[ctl,:] = el
- ctl+=1
- if Bords['Left']=='Abso':
- Elements2DBordAbso[cta,:] = el
- cta+=1
- if Bords['Left']=='Free':
- Elements2DBordFree[ctf,:] = el
- ctf+=1
- # Cas du bord double
- if len(rr)==3:
- ctdb=0
- # Recherche des deux elements 1D
- for Nn in xrange(Ninc1DBordLeft):
- kk = set(Elements1DBordLeft[Nn])
- if not set.isdisjoint(rr, kk):
- if kk.issubset(jj):
- exec 'kk'+str(ctdb)+'=kk'
- ctdb+=1
- Elc = set.intersection(kk0, kk1)
- rrm=rr.difference(Elc)
- rr1=list(rrm)
- El1=concatenate(([rr1[0]],list(Elc)))
- El2=concatenate((list(Elc),[rr1[1]]))
- el1 = concatenate(([Ct2D+1], [2], El1))
- el2 = concatenate(([Ct2D+1], [2], El2))
- Elements2DBordLeft[ctl,:] = el1
- ctl+=1
- Elements2DBordLeft[ctl,:] = el2
- ctl+=1
- if Bords['Left']=='Abso':
- Elements2DBordAbso[cta,:] = el1
- cta+=1
- Elements2DBordAbso[cta,:] = el2
- cta+=1
- if Bords['Left']=='Free':
- Elements2DBordFree[ctf,:] = el1
- ctf+=1
- Elements2DBordFree[ctf,:] = el2
- ctf+=1
- # Bord Top
- if not set.isdisjoint(jj, NodesBordTop2n):
- rr = set.intersection(jj, NodesBordTop2n)
- if len(rr)==2:
- el = concatenate(([Ct2D+1], [len(rr)], list(rr)))
- Elements2DBordTop[ctt,:len(el)] = el
- ctt+=1
- if Bords['Top']=='Abso':
- Elements2DBordAbso[cta,:] = el
- cta+=1
- if Bords['Top']=='Free':
- Elements2DBordFree[ctf,:] = el
- ctf+=1
- # Bord Bottom
- if not set.isdisjoint(jj, NodesBordBottom2n):
- rr = set.intersection(jj, NodesBordBottom2n)
- if len(rr)==2:
- el = concatenate(([Ct2D+1], [len(rr)], list(rr)))
- Elements2DBordBottom[ctb,:len(el)] = el
- ctb+=1
- if Bords['Bottom']=='Abso':
- Elements2DBordAbso[cta,:] = el
- cta+=1
- if Bords['Bottom']=='Free':
- Elements2DBordFree[ctf,:] = el
- ctf+=1
- Elements2DBord=Elements2DBord[:ctBord]
+ print '-'*60
+ print 'Number of elements in contact with a free surface', ctf
+ print 'Number of elements in contact with an absorbing surface', cta
+ print '#'*20+' END OF PROCESSING '+'#'*20
+ # remove unnecessary zeros created at initialization
#----------------------------------------------------------------------
Elements2DBordAbso = Elements2DBordAbso[:cta,:]
Elements2DBordFree = Elements2DBordFree[:ctf,:]
- Elements2DBordTop = Elements2DBordTop[:ctt,:]
- Elements2DBordBottom = Elements2DBordBottom[:ctb,:]
- Elements2DBordRight = Elements2DBordRight[:ctr,:]
- Elements2DBordLeft = Elements2DBordLeft[:ctl,:]
#-----------------------------------------------------------------------
- # Sauvegarde au format SPECFEM
- SauvFicSpecfem('Mesh_'+Nom, Ninc2D, Elements, '%i')
+ # Save to SPECFEM format
+ SauvFicSpecfem('Mesh_'+Nom, N2D, Elements, '%i')
#
SauvFicSpecfem('Surf_abs_'+Nom, cta, Elements2DBordAbso, '%i')
#
@@ -368,22 +377,20 @@
#
savetxt('Material_'+Nom,Milieu, fmt='%i')
#
- SauvFicSpecfem('Surf_top_'+Nom, ctt, Elements2DBordTop, '%i')
+ if PML: SauvFicSpecfem('EltPML_'+Nom, N2DPML, ElementsPML, '%i')
#
- SauvFicSpecfem('Surf_bottom_'+Nom, ctb, Elements2DBordBottom, '%i')
- #
- SauvFicSpecfem('Surf_right_'+Nom, ctr, Elements2DBordRight, '%i')
- #
- SauvFicSpecfem('Surf_left_'+Nom, ctl, Elements2DBordLeft, '%i')
return
if __name__=='__main__':
- set_printoptions(precision=6, threshold=None, edgeitems=None, linewidth=200, suppress=None, nanstr=None, infstr=None)
+ set_printoptions(precision=6, threshold=None, edgeitems=None, \
+ linewidth=200, suppress=None, nanstr=None, infstr=None)
#
- # Lecture des paramètres d'entree eventuels
+ # Reading input parameter if provided
#
Fic = sys.argv[1]; del sys.argv[1]
#
- Bords={'Top':'Abso', 'Bottom':'Abso', 'Left':'Abso' , 'Right':'Abso' }
+ # Default values
+ Bords = {'Top':'Abso', 'Bottom':'Abso', 'Left':'Abso' , 'Right':'Abso' }
+ #--------------------------------------------------------------------------
while len(sys.argv) > 1:
opt=sys.argv[1]; del sys.argv[1]
if opt == '-t':
@@ -410,7 +417,7 @@
else:
print 'Wrong condition'
del sys.argv[1]
- elif opt == '-r': # On affiche le resultat ou pas
+ elif opt == '-r':
if sys.argv[1]== 'F':
Bords['Right']='Free'
elif sys.argv[1]== 'A':
More information about the CIG-COMMITS
mailing list