forked from asaintenoy/Porchet-GPR
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ecriture_fichiers_GPRMAX.py
123 lines (99 loc) · 5.98 KB
/
ecriture_fichiers_GPRMAX.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
import numpy as np
import math
def ecriture_fichiers_GPRMAX(X, Y, grid_z0, trace_number, nom, paramMVG, paramGPRMAX, geometry, dl, materiaux) :
#attention, veut les param en m!!!!
etrou = geometry.etrou*0.01
h_eau = geometry.h_eau*0.01
radius = geometry.r*0.01
Ks = paramMVG.Ks
sigma = paramGPRMAX.sigma
nmedia=len(materiaux)+2
dx=abs((X)[0,1]-(X)[0,0])
dy=abs((Y)[1,0]-(Y)[0,0])
fgrid=open(nom+'_'+str(trace_number+1)+'.in',"w")
fgrid.write("""------------------------------------------------\n""")
#fgrid.write("""#number_of_media: {}\n""".format(nmedia+4))
fgrid.write("""#domain: {} {} {}\n""".format(round(2*max(X[0,:])+dx,2), round(max(Y[:,0])*2,2), dl))
fgrid.write("""#dx_dy_dz: {} {} {}\n""".format(dl, dl, dl))
fgrid.write("""#time_window: {}\n""".format(paramGPRMAX.time))
#fgrid.write("""#time_step_stability_factor: {}\n""".format(fac_dt))
#fgrid.write("""#abc_type: pml\n""")
#fgrid.write("""#pml_cells: {}\n""".format(10))
fgrid.write("""#waveform: ricker 1.0 {} Mydipole\n""".format(paramGPRMAX.wave_freq))
w=max(X[0,:])+dx/2 #largeur
fgrid.write("""#hertzian_dipole: z {} {} {} Mydipole\n""".format(round(w+paramGPRMAX.d_emet,2), round(max(Y[:,0])+3/2*dy,2), 0.0))#TX
fgrid.write("""#rx: {} {} {} \n""".format(round(w+paramGPRMAX.d_recept,2), round(max(Y[:,0]+3/2*dy),2), 0.0))#RX
fgrid.write("""------------------------------------------------\n""")
k=1
VWC=np.zeros([len(Y[:,0])*len(grid_z0[0,:])+1, 1])
for i in range(0, len(Y[:,0])) :
for j in range(0, len(grid_z0[0,:])) :
VWC[k-1]=grid_z0[i,j]
k = k+1
#SI EAU DOUCE INJECTEE
for i in materiaux:
fgrid.write("""#material: {} {} 1.0 0.0 {}\n""".format(i, sigma, materiaux[i][0]))
#SI EAU SALEE INJECTEE
#TODO: Mettre un test pour demander à l'utilisateur quel cas il veut faire...
#for i in materiaux:
# fgrid.write("""#material: {} {} 1.0 0.0 {}\n""".format(i, round(materiaux[i][1],2), materiaux[i][0]))
#SI MATERIAUX EAU ET/OU PVC PRESENTS DANS LE MODELE
#fgrid.write("""#material: {} {} 1.0 0.0 eau\n""".format(paramGPRMAX.eps_w,sigma)) #Donne le epsilon de l'eau
#fgrid.write("""#material: {} {} 1.0 0.0 pvc\n""".format(paramGPRMAX.eps_pvc,sigma)) #Donne le epsilon du pvc
fgrid.write("""------------------------------------------------\n""")
#Creation du modele en entier par symetrie des resultats SWMS
print(trace_number+1)
d=dy/2 #hauteur
k=1 #compteur
count=0
A=np.zeros([len(X[:,0])*len(X[0,:])*2, 3]) #Taille de la matrice = nombre de cellules du nouveau maillage*2
A_tab={}
for ii in range(0,len(Y[:,0])) :
for jj in range(0, len(grid_z0[0,:])) :
fgrid.write("""#box: {} {} {} {} {} {} {}\n""".format(round(w+X[ii,jj]-dx/2,2), round(d+Y[ii,jj]-dy/2,2), 0.0, round(w+X[ii,jj]+dx/2,2),round(d+Y[ii,jj]+dy/2,2), dl, materiaux[grid_z0[ii,jj]][0]))
A[count,0]=w+X[ii,jj]-dx/2
A[count,1]=d+Y[ii,jj]-dy/2
A[count,2]=grid_z0[ii,jj]
count=count+1
fgrid.write("""#box: {} {} {} {} {} {} {}\n""".format(round(w-X[ii,jj]-dx/2,2), round(d+Y[ii,jj]-dy/2,2), 0.0, round(w-X[ii,jj]+dx/2,2),round(d+Y[ii,jj]+dy/2,2), dl, materiaux[grid_z0[ii,jj]][0]))
A[count,0]=w-X[ii,jj]-dx/2
A[count,1]=d+Y[ii,jj]-dy/2
A[count,2]=grid_z0[ii,jj]
k=k+1
count=count+1
#On bouche le tuyau et on interpole la fermeture du bulbe
fgrid.write("""------------------------------------------------\n""")
o=np.where((A[:,0]<=w+dx/2+radius) & (A[:,0]>=w-dx/2-radius) & (A[:,1]>=etrou))
ii=np.unique(A[o][:,1])
blou=A[o]
for i in ii:
a=np.where(A[o,1]==i)
b=np.where((A[:,0]==w+dx/2+radius) & (A[:,1]==i))
for ii in a[1]:
blou[ii,2]=A[b,2]
A[o]=blou
fgrid.write("""#box: {} {} {} {} {} {} {}\n""".format(round(blou[ii,0],2), round(i,2),0.0, round(blou[ii,0]+dx,2), round(i+dy,2), dl, materiaux[blou[ii,2]][0]))
#Si premier pas de temps, on ajoute une boite d'eau
if trace_number==0 :
q=np.where((A[:,0]<=w+dx/2+radius) & (A[:,0]>=w-dx/2-radius) & (A[:,1]>=etrou) & (A[:,1]<=etrou+h_eau))
A[q,2]=grid_z0.max()
fgrid.write("""#box: {} {} {} {} {} {} {}\n""".format(w-dx/2-radius, etrou, 0.0, w+dx/2+radius, etrou+h_eau, dl, materiaux[grid_z0.max()][0])) #materiaux[grid_z0.max()][0]
#Ajout du tuyau pvc dans le trou
#fgrid.write("""#box: {} {} {} {} {} {} pvc\n""".format(round(w-dx/2-radius,2), round(etrou+d,2), 0.0, round(w+dx/2+radius,2), round(max(Y[:,0])*1.1,2), dl)) #ajout du tuyau tout le long
#fgrid.write("""#box: {} {} {} {} {} {} pvc\n""".format(0.38, 0.505, 0.0, 0.43, 0.88, 0.6*dl)) #ajout du tuyau tout le long
#o=np.where((A[:,0]<=w+dx/2+radius) & (A[:,0]>=w-dx/2-radius) & (A[:,1]>=etrou))
#o=np.where((A[:,0]<=w+dx/2+radius) & (A[:,0]>=w-dx/2-radius) & (A[:,1]>=etrou))
#A[o,2]=eps_pvc
#Ajout d'eau au fond du trou
#fgrid.write("""#box: {} {} {} {} {} {} eau\n""".format(w-dx/2-radius+0.007, etrou+d, 0.0, w+dx/2+radius-0.007, etrou+d+h_eau, dl))
#q=np.where((A[:,0]<=w+dx/2+radius-0.007) & (A[:,0]>=w-dx/2-radius+0.007) & (A[:,1]<=etrou+d+h_eau) & (A[:,1]>=etrou+d))
#A[q,2]=81
#Ajout de l'air dans le trou
#fgrid.write("""#box: {} {} {} {} {} {} free_space\n""".format(w-dx/2-radius+0.007, etrou, 0.0, w+dx/2+radius-0.007, max(Y[:,0])*2, dl))
#s=np.where((A[:,0]<=w+dx/2+radius-0.007) & (A[:,0]>=w-dx/2-radius+0.007) & (A[:,1]<=max(Y[:,0])*1.1) & (A[:,1]>=etrou+d+h_eau))
#A[s,2]=1
fgrid.write("""------------------------------------------------\n""")
fgrid.write("""#messages: n\n""")
#fgrid.write("""#geometry_view: 0 0 0 {} {} {} {} {} {} modele_vue n""".format(2*max(X[0,:])+dx, max(Y[:,0])*2, dl, dl, dl, dl))
fgrid.close()
return A