Commit 48b7be4b authored by Kévin Moreau's avatar Kévin Moreau

Add scripts to run several test cases at once

parent dca1323f
#!/usr/bin/python2
import matplotlib.pyplot as plt
import numpy as np
# dirs
loc="/glouton/struct/kmoreau/exlibris/test/relduct_dam_hole/"
# data
radii=np.linspace(0., 0.75, 16)
lcs=np.linspace(0.15, 0.55, 5)
for lc in lcs:
Fmaxs=[]
F0s=[]
plt.figure(1)
plt.xlabel('U/U_c')
plt.ylabel('F/F_c')
for r in radii:
jobname="radius2_"+str(int(100*r))+"_lc_"+str(int(100*lc))
resloc=loc+jobname+"/"
in1=np.loadtxt(resloc+"sensor_val_disp_max.txt", delimiter=' ', skiprows=0)
in3=np.loadtxt(resloc+"sensor_val_force_y.txt", delimiter=' ', skiprows=0)
in2=np.array([[0.,0.],[in1[0,1],in3[0,1]]])
# plt.plot(in2[:,0]/in1[0,1], in2[:,1]/in3[0,1], '+-')
# plt.plot(in1[:,1]/in1[0,1], in3[:,1]/in3[0,1], '+-')
Fmaxs.append(np.max(in3[:,1]))
F0s.append(in3[0,1])
plt.figure(2)
plt.plot(radii, Fmaxs/Fmaxs[0], '*-', label="lc="+str(lc))
plt.plot(radii, 1./3.*np.ones(len(radii)), '--')
plt.xlabel('r')
plt.ylabel('Fmax/Fmax_0')
plt.legend(loc=1)
plt.show()
#!/usr/bin/python2
import re
import subprocess as sp
import numpy as np
# dirs
dataloc="/home/users/struct/kmoreau/code/exlibris/Applis/TLSDuctile/test/dam_hole/data/"
loc="/glouton/struct/kmoreau/exlibris/test/relduct_dam_hole/"
# dataloc="/home/kev/documents/code/exlibris/Applis/TLSDuctile/test/dam_hole/data/"
# loc="/scratch/exlibris/test/relduct_dam_hole/"
# data
radii=np.linspace(0., 0.75, 16)
lcs=np.linspace(0.15, 0.55, 5)
f=open(dataloc+"radii_template.geo", 'r')
geodata=f.read()
f.close()
f=open(dataloc+"info.dat", 'r')
infodata=f.read()
f.close()
for r in radii:
for lc in lcs:
jobname="radius2_"+str(int(100*r))+"_lc_"+str(int(100*lc))
resloc=loc+jobname+"/"
# create result dir
sp.call(["mkdir", resloc])
# import data files
sp.call(["cp", dataloc+"main.dat", resloc+"main.dat"])
sp.call(["cp", dataloc+"elasto_dam.mat", resloc+"mate.mat"])
# correct data files
infodata=re.sub("lc = [0-9.]+", "lc = "+str(lc), infodata)
filename=resloc+"info.dat"
f=open(filename, 'w')
f.write(infodata)
f.close()
# create mesh
geodata=re.sub("r=[0-9.e\-]+", "r="+str(r), geodata)
filename=resloc+"mesh.geo"
f=open(filename, 'w')
f.write(geodata)
f.close()
sp.call(["gmsh", "-2", "-format", "msh1", filename])
# create script
filename=resloc+jobname+".sh"
f=open(filename, 'w')
f.write("#!/bin/sh\n")
f.write("cd "+resloc+"\n")
f.write("time ../dam_hole main.dat &> out.log\n")
f.close()
sp.call(["chmod", "+x", filename])
sp.call(["qsub", "-qdefaut", "-lncpus=1,mem=4gb", filename])
#!/usr/bin/python2
import matplotlib.pyplot as plt
import numpy as np
# dirs
loc="/glouton/struct/kmoreau/exlibris/test/relduct_dam_hole/"
# data
radii=np.linspace(0., 0.36, 13)
lcs=np.linspace(0.15, 0.55, 5)
for lc in lcs:
Fmaxs=[]
F0s=[]
plt.figure(1)
plt.xlabel('U/U_c')
plt.ylabel('F/F_c')
for r in radii:
jobname="radius_"+str(int(100*r))+"_lc_"+str(int(100*lc))
resloc=loc+jobname+"/"
in1=np.loadtxt(resloc+"sensor_val_disp_max.txt", delimiter=' ', skiprows=0)
in3=np.loadtxt(resloc+"sensor_val_force_y.txt", delimiter=' ', skiprows=0)
in2=np.array([[0.,0.],[in1[0,1],in3[0,1]]])
# plt.plot(in2[:,0]/in1[0,1], in2[:,1]/in3[0,1], '+-')
# plt.plot(in1[:,1]/in1[0,1], in3[:,1]/in3[0,1], '+-')
Fmaxs.append(np.max(in3[:,1]))
F0s.append(in3[0,1])
plt.figure(2)
plt.plot(radii, Fmaxs/Fmaxs[0], '*-', label="lc="+str(lc))
plt.plot(radii, 1./3.*np.ones(len(radii)), '--')
plt.xlabel('r')
plt.ylabel('Fmax/Fmax_0')
plt.legend(loc=1)
plt.show()
#!/usr/bin/python2
import re
import subprocess as sp
import numpy as np
# dirs
dataloc="/home/users/struct/kmoreau/code/exlibris/Applis/TLSDuctile/test/dam_hole/data/"
loc="/glouton/struct/kmoreau/exlibris/test/relduct_dam_hole/"
# dataloc="/home/kev/documents/code/exlibris/Applis/TLSDuctile/test/dam_hole/data/"
# loc="/scratch/exlibris/test/relduct_dam_hole/"
# data
radii=np.linspace(0., 0.36, 13)
lcs=np.linspace(0.15, 0.55, 5)
f=open(dataloc+"radii_template.geo", 'r')
geodata=f.read()
f.close()
f=open(dataloc+"info.dat", 'r')
infodata=f.read()
f.close()
for r in radii:
for lc in lcs:
jobname="radius_"+str(int(100*r))+"_lc_"+str(int(100*lc))
resloc=loc+jobname+"/"
# create result dir
sp.call(["mkdir", resloc])
# import data files
sp.call(["cp", dataloc+"main.dat", resloc+"main.dat"])
sp.call(["cp", dataloc+"elasto_dam.mat", resloc+"mate.mat"])
# correct data files
infodata=re.sub("lc = [0-9.]+", "lc = "+str(lc), infodata)
filename=resloc+"info.dat"
f=open(filename, 'w')
f.write(infodata)
f.close()
# create mesh
geodata=re.sub("r=[0-9.e\-]+", "r="+str(r), geodata)
filename=resloc+"mesh.geo"
f=open(filename, 'w')
f.write(geodata)
f.close()
sp.call(["gmsh", "-2", "-format", "msh1", filename])
# create script
filename=resloc+jobname+".sh"
f=open(filename, 'w')
f.write("#!/bin/sh\n")
f.write("cd "+resloc+"\n")
f.write("time ../dam_hole main.dat &> out.log\n")
f.close()
sp.call(["chmod", "+x", filename])
sp.call(["qsub", "-qdefaut", "-lncpus=1,mem=2gb", filename])
#!/usr/bin/python2
import matplotlib.pyplot as plt
import numpy as np
N=5000
eps=np.zeros(N)
sig=np.zeros(N)
d=np.zeros(N)
# hardening function and derivative
def h(d):
return 1./np.power(1.-d, 1./n)
def dh(d):
return 1./n/np.power(1.-d, (n+1.)/n)
def hrec(e):
return 1.-np.power(1./e, n)
def compute_local_law_explicit():
eps[0]=0
sig[0]=0
d[0]=0
for i in range(1, N):
eps[i]=eps[i-1]+deps
if eps[i]>eps_c:
dd=E*eps[i]*deps/(Y_c*dh(d[i-1]))
d[i]=min(d[i-1]+dd,1)
sig[i]=max(0,(1-d[i])*E*eps[i])
# main
E=210.e9
Y_c=1.32e5
eps_c=np.sqrt(2*Y_c/E)
sig_c=E*eps_c
ns=[1]
deps=eps_c/1000.
for n in ns:
compute_local_law_explicit()
plt.plot(eps/eps_c, sig/sig_c, '--')
m1=eps<1.001*eps_c
m2=eps>eps_c
plt.plot(eps[m1]/eps_c, E*eps[m1]/sig_c)
plt.plot(eps[m2]/eps_c, (1-hrec(np.power(eps[m2]/eps_c, 2.)))*E*eps[m2]/sig_c)
# for i in ["run"]:
# in1=np.loadtxt(i+"/sensor_val_disp_max.txt", delimiter=' ', skiprows=0)
# in3=np.loadtxt(i+"/sensor_val_force_y.txt", delimiter=' ', skiprows=0)
# in2=np.array([[0.,0.],[in1[0,1],in3[0,1]]])
# plt.plot(in2[:,0]/in1[0,1], in2[:,1]/in3[0,1], '+')
# plt.plot(in1[:,1]/in1[0,1], in3[:,1]/in3[0,1], '+')
plt.show()
import matplotlib.pyplot as plt
import numpy as np
for i in ["dupl02", "dupl07"]:
in1=np.loadtxt(i+"/sensor_val_disp_max.txt", delimiter=' ', skiprows=0)
in3=np.loadtxt(i+"/sensor_val_force_y.txt", delimiter=' ', skiprows=0)
in2=np.array([[0.,0.],[in1[0,1],in3[0,1]]])
coeff=in1[0,1]
plt.plot(in2[:,0]/coeff, in2[:,1]/in3[0,1], '+-')
plt.plot(in1[:,1]/coeff, in3[:,1]/in3[0,1], '+-')
plt.xlim(0, 4)
plt.ylabel('T')
plt.xlabel('U')
plt.show()
#!/usr/bin/python2
import matplotlib.pyplot as plt
import numpy as np
# dirs
loc="/glouton/struct/kmoreau/exlibris/test/relduct_dam_notch/"
# data
angles=np.linspace(0, 90, 19)
angles[0]=1
lcs=np.linspace(0.15, 0.55, 5)
for lc in lcs:
Fmaxs=[]
F0s=[]
plt.figure(1)
plt.xlabel('U/U_c')
plt.ylabel('F/F_c')
for theta in angles:
jobname="angle2_"+str(int(theta))+"_lc_"+str(int(100*lc))
resloc=loc+jobname+"/"
in1=np.loadtxt(resloc+"sensor_val_disp_max.txt", delimiter=' ', skiprows=0)
in3=np.loadtxt(resloc+"sensor_val_force_y.txt", delimiter=' ', skiprows=0)
in2=np.array([[0.,0.],[in1[0,1],in3[0,1]]])
# plt.plot(in2[:,0]/in1[0,1], in2[:,1]/in3[0,1], '+-')
# plt.plot(in1[:,1]/in1[0,1], in3[:,1]/in3[0,1], '+-')
Fmaxs.append(np.max(in3[:,1]))
F0s.append(in3[0,1])
plt.figure(2)
plt.plot(angles, Fmaxs, '*-', label="lc="+str(lc))
plt.xlabel('theta')
plt.ylabel('Fmax')
plt.legend(loc=2)
plt.show()
#!/usr/bin/python2
import re
import subprocess as sp
import numpy as np
# dirs
dataloc="/home/users/struct/kmoreau/code/exlibris/Applis/TLSDuctile/test/dam_notch/data/"
loc="/glouton/struct/kmoreau/exlibris/test/relduct_dam_notch/"
# dataloc="/home/kev/documents/code/exlibris/Applis/TLSDuctile/test/dam_notch/data/"
# loc="/scratch/exlibris/test/relduct_dam_notch/"
# data
angles=np.linspace(0, 90, 19)
angles[0]=1
lcs=np.linspace(0.15, 0.55, 5)
Gc=40000
f=open(dataloc+"angles_template.geo", 'r')
geodata=f.read()
f.close()
f=open(dataloc+"info.dat", 'r')
infodata=f.read()
f.close()
f=open(dataloc+"elasto_dam.mat", 'r')
matedata=f.read()
f.close()
for theta in angles:
for lc in lcs:
jobname="angle2_"+str(int(theta))+"_lc_"+str(int(100*lc))
resloc=loc+jobname+"/"
# create result dir
sp.call(["mkdir", resloc])
# import data files
sp.call(["cp", dataloc+"main.dat", resloc+"main.dat"])
# correct data files
infodata=re.sub("lc = [0-9.]+", "lc = "+str(lc), infodata)
filename=resloc+"info.dat"
f=open(filename, 'w')
f.write(infodata)
f.close()
matedata=re.sub("Y_CRIT = [0-9.e\-]+", "Y_CRIT = "+str(Gc/lc), matedata)
filename=resloc+"mate.mat"
f=open(filename, 'w')
f.write(matedata)
f.close()
# create mesh
radtheta=theta*(2.*np.pi/360.)
sintheta=np.sin(radtheta)
costheta=np.cos(radtheta)
geodata=re.sub("sintheta=[0-9.e\-]+", "sintheta="+str(sintheta), geodata)
geodata=re.sub("costheta=[0-9.e\-]+", "costheta="+str(costheta), geodata)
filename=resloc+"mesh.geo"
f=open(filename, 'w')
f.write(geodata)
f.close()
sp.call(["gmsh", "-2", "-format", "msh1", filename])
# create script
filename=resloc+jobname+".sh"
f=open(filename, 'w')
f.write("#!/bin/sh\n")
f.write("cd "+resloc+"\n")
f.write("time ../dam_notch main.dat &> out.log\n")
f.close()
sp.call(["chmod", "+x", filename])
sp.call(["qsub", "-qdefaut", "-lncpus=1,mem=4gb", filename])
#!/usr/bin/python2
import matplotlib.pyplot as plt
import numpy as np
# dirs
loc="/glouton/struct/kmoreau/exlibris/test/relduct_dam_notch/"
# data
angles=np.linspace(0, 90, 19)
angles[0]=1
# angles=[45]
ns=np.array([2, 3, 4, 5, 6, 7, 8])
# ns=np.array([2])
for n in ns:
Fmaxs=[]
F0s=[]
plt.figure(1)
plt.xlabel('U/U_c')
plt.ylabel('F/F_c')
for theta in angles:
jobname="angle3_2nd_"+str(int(theta))+"_n_"+str(int(n))
resloc=loc+jobname+"/"
in1=np.loadtxt(resloc+"sensor_val_disp_max.txt", delimiter=' ', skiprows=0)
in3=np.loadtxt(resloc+"sensor_val_force_y.txt", delimiter=' ', skiprows=0)
in2=np.array([[0.,0.],[in1[0,1],in3[0,1]]])
# plt.plot(in2[:,0]/in1[0,1], in2[:,1]/in3[0,1], '+-')
# plt.plot(in1[:,1]/in1[0,1], in3[:,1]/in3[0,1], '+-')
Fmaxs.append(np.max(in3[:,1]))
F0s.append(in3[0,1])
plt.figure(2)
plt.plot(angles, Fmaxs, '*-', label="n="+str(n))
plt.xlabel('theta')
plt.ylabel('Fmax')
plt.legend(loc=2)
plt.show()
#!/usr/bin/python2
import re
import subprocess as sp
import numpy as np
# dirs
dataloc="/home/users/struct/kmoreau/code/exlibris/Applis/TLSDuctile/test/dam_notch/data/"
loc="/glouton/struct/kmoreau/exlibris/test/relduct_dam_notch/"
# dataloc="/home/kev/documents/code/exlibris/Applis/TLSDuctile/test/dam_notch/data/"
# loc="/scratch/exlibris/test/relduct_dam_notch/"
# data
angles=np.linspace(0, 90, 19)
angles[0]=1
ns=np.array([2, 3, 4, 5, 6, 7, 8])
f=open(dataloc+"angles_template.geo", 'r')
geodata=f.read()
f.close()
f=open(dataloc+"info.dat", 'r')
infodata=f.read()
f.close()
infodata=re.sub("lc = 0.5", "lc = 0.35", infodata)
f=open(dataloc+"elasto_dam.mat", 'r')
matedata=f.read()
f.close()
matedata=re.sub("HARDENING = exponential", "HARDENING = rational", matedata)
matedata=re.sub("HARDENING_COEFF = 4.", "HARDENING_COEFF = 1.", matedata)
m=re.search("Y_CRIT = [0-9.e\-]+", matedata)
Yc=m.group(0).split()[2]
for theta in angles:
for n in ns:
jobname="angle3_"+str(int(theta))+"_n_"+str(int(n))
resloc=loc+jobname+"/"
# create result dir
sp.call(["mkdir", resloc])
# import data files
sp.call(["cp", dataloc+"main.dat", resloc+"main.dat"])
# correct data files
filename=resloc+"info.dat"
f=open(filename, 'w')
f.write(infodata)
f.close()
matedata=re.sub("Y_CRIT = [0-9.e\-]+", "Y_CRIT = "+str(Yc*(1.-1./n)), matedata)
filename=resloc+"mate.mat"
f=open(filename, 'w')
f.write(matedata+"HARDENING_POWER = "+str(1./n)+"\n")
f.close()
# create mesh
radtheta=theta*(2.*np.pi/360.)
sintheta=np.sin(radtheta)
costheta=np.cos(radtheta)
geodata=re.sub("sintheta=[0-9.e\-]+", "sintheta="+str(sintheta), geodata)
geodata=re.sub("costheta=[0-9.e\-]+", "costheta="+str(costheta), geodata)
filename=resloc+"mesh.geo"
f=open(filename, 'w')
f.write(geodata)
f.close()
sp.call(["gmsh", "-2", "-format", "msh1", filename])
# create script
filename=resloc+jobname+".sh"
f=open(filename, 'w')
f.write("#!/bin/sh\n")
f.write("cd "+resloc+"\n")
f.write("time ../dam_notch main.dat &> out.log\n")
f.close()
sp.call(["chmod", "+x", filename])
sp.call(["qsub", "-qdefaut", "-lncpus=1,mem=4gb", filename])
#!/usr/bin/python2
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter
# dirs
loc='/glouton/struct/kmoreau/exlibris/test/relduct_dam_notch/'
suffix='_last'
# suffix=''
# data
angles=np.linspace(0., 90., 19)
angles[0]=1.
angles[18]=89.
N=10
DeltaX=(np.log(300.)-np.log(3.))/N
X=np.log(3.)+np.linspace(0., 10., N+1)*DeltaX
heights=np.around(np.exp(X), decimals=2)
angles_curve=[45]
angles_log=np.linspace(0., 90., 10)
angles_log[0]=1.
angles_log[9]=89.
i=1
for theta in angles_curve:
for h in heights:
jobname='angle_'+str(int(theta))+'_height_'+str(int(h))+suffix
resloc=loc+jobname+'/'
in1=np.loadtxt(resloc+'sensor_val_load_factor.txt')
in3=np.loadtxt(resloc+'sensor_val_force_y.txt')
in2=np.array([[0.,0.],[in1[0,1],in3[0,1]]])
coeff=h/100.
# coeff=1.
plt.figure(i)
plt.plot(in2[:,0]*coeff, in2[:,1], '-')
plt.plot(in1[:,1]*coeff, in3[:,1], '-', label='$h/l_c='+str(h)+'$')
plt.figure(i+1)
plt.plot(in2[:,0]*coeff/h, in2[:,1]/h, '-')
plt.plot(in1[:,1]*coeff/h, in3[:,1]/h, '-', label='$h/l_c='+str(h)+'$')
if theta==90:
jobname='angle_90_height_infty'+suffix
resloc=loc+jobname+'/'
in1=np.loadtxt(resloc+'sensor_val_load_factor.txt')
in3=np.loadtxt(resloc+'sensor_val_force_y.txt')
in2=np.array([[0.,0.],[in1[0,1],in3[0,1]]])
plt.figure(i+1)
plt.plot(in2[:,0]/h, in2[:,1]/h, '-')
plt.plot(in1[:,1]/h, in3[:,1]/h, '-', label='$h/l_c=+\infty$')
i+=2
for h in heights:
force_maxs=[]
jobname='angle_89_height_'+str(int(h))+suffix
resloc=loc+jobname+'/'
in3=np.loadtxt(resloc+'sensor_val_force_y.txt')
force_y_max_90=np.amax(in3[:,1])
for theta in angles:
jobname='angle_'+str(int(theta))+'_height_'+str(int(h))+suffix
resloc=loc+jobname+'/'
in3=np.loadtxt(resloc+'sensor_val_force_y.txt')
force_maxs.append(np.amax(in3[:,1]))
plt.figure(i)
plt.plot(angles, np.array(force_maxs)/force_y_max_90, '+-', label='$h/l_c='+str(h)+'$')
i+=1
# log(L/lc) vs log(sig_max^theta/sig_max^90)
for theta in angles_log:
absc=[]
ordo=[]
for h in heights:
jobname='angle_89_height_'+str(int(h))+suffix
resloc=loc+jobname+'/'
in3=np.loadtxt(resloc+'sensor_val_force_y.txt')
force_y_max_90=np.amax(in3[:,1])
jobname='angle_'+str(int(theta))+'_height_'+str(int(h))+suffix
resloc=loc+jobname+'/'
in3=np.loadtxt(resloc+'sensor_val_force_y.txt')
absc.append(h)
ordo.append(np.amax(in3[:,1]/force_y_max_90))
plt.figure(i)
plt.plot(absc, ordo, '+-', label='$\\theta='+str(theta)+'$')
# plotting
i=1
for theta in angles_curve:
plt.figure(i)
plt.xlabel('$U$')
plt.ylabel('$F$')
plt.xlim([0., 0.075])
plt.legend(loc=1)
plt.figure(i+1)
plt.xlabel('$U/h$')
plt.ylabel('$F/h$')
plt.xlim([0., 3.e-3])
plt.legend(loc=1)
i+=2
plt.figure(i)
plt.xlabel('$\\theta$')
plt.ylabel('$\sigma_{\mathsf{max}}^{\\theta}/\sigma_{\mathsf{max}}^{89}$')
plt.ylim([0.1, 1.1])
plt.legend(loc=1)
i+=1
fig=plt.figure(i)
ax=fig.add_subplot(1,1,1)
x=np.linspace(3., 300., 30)
X=304
Y=0.120
a=-1./2.
plt.plot(x, Y*np.power(x/X, a), color='g', linestyle='--', label="-1/2")
# plt.xlim([3.5, 13.5])
plt.ylim([0.1, 1.1])
plt.xscale('log', basex=np.e)
plt.yscale('log', basey=np.e)
plt.xlabel('$h/l_c$')
plt.ylabel('$\sigma_\mathsf{max}^\\theta/\sigma_\mathsf{max}^{89}$')
plt.legend(loc=3)
ax.xaxis.set_major_formatter(ScalarFormatter())
ax.yaxis.set_major_formatter(ScalarFormatter())
for fig in range(1, i+1):
plt.figure(fig)
plt.legend(loc=2, bbox_to_anchor=(1.05,1), ncol=2)
plt.savefig('fig'+str(fig)+'.pdf')
# plt.show()
#!/usr/bin/python2
import re
import subprocess as sp
import numpy as np
import time
import math
# dirs
# dataloc="/home/users/struct/kmoreau/code/exlibris/Applis/TLSDuctile/test/dam_notch/data/"
# loc="/glouton/struct/kmoreau/exlibris/test/relduct_dam_notch/"
dataloc="/home/kev/documents/code/exlibris/Applis/TLSDuctile/test/dam_notch/data/"
loc="/scratch/exlibris/test/relduct_dam_notch/"
suffix='_adou'
# data
angles=np.linspace(0., 90., 19)
angles[0]=1.
angles[18]=89.
N=10
DeltaX=(np.log(600.)-np.log(6.))/N