Commit a68a8e41 authored by Gilles MARCKMANN's avatar Gilles MARCKMANN

Merge branch 'master' of git.gem.ec-nantes.fr:TTK/TLSDuctile

parents af8e97de 75e5050b
Pipeline #243 skipped
......@@ -206,8 +206,9 @@ double QSFormulation::computeLoadFactor() {
std::cout<<"load factor "<<load_factor<<" ref force norm "<<ref_force_norm<< " force norm "<<load_factor*ref_force_norm<<" ref disp norm "<<ref_disp_norm<<" disp norm "<<load_factor*ref_disp_norm<<std::endl;
if(isPeakForce()) {
OldAndCurrent_c::old();
post_pro.exportOnSpace("damage_peak", 0, tls_solver.getDamageEval(), geom.getIntegRulePartition(0), geom.begin(), geom.end());
post_pro.exportOnSpace("damage_peak", 0, tls_solver.getDamageEval(), geom.getIntegRuleSmart(), geom.begin(), geom.end());
OldAndCurrent_c::current();
post_pro.exportOnTime("peak_step", Observer::tell<int>("step"), 1.);
}
return load_factor;
}
......
......@@ -132,8 +132,7 @@ double ExponentialHardeningFunction::getGrad(double x) const {
return c*exp(c*x);
}
double ExponentialHardeningFunction::getPrimitive(double x) const {
std::cout<<"TODO"<<std::endl;
return 0.;
return exp(c*x)/c;
}
RationalHardeningFunction::RationalHardeningFunction(const std::string& filename) {
......
......@@ -385,6 +385,10 @@ void TLSGeom::computeNonlocalBoundingBox() const {
post_pro.exportOnTime("nonlocal_bb_max_x", step, max(0));
post_pro.exportOnTime("nonlocal_bb_max_y", step, max(1));
post_pro.exportOnTime("nonlocal_bb_max_z", step, max(2));
// if(max(0)>4.e-4) {
// std::abort();
// }
}
void TLSGeom::updateNonlocal(const std::list<mEntity*>& delocalized_list) {
......
......@@ -274,8 +274,6 @@ void TLSSolver::delocalizeLevelSetField() {
LocateHotPointsCommand locate_hot_points(eval_norm_grad_phi, 1., delocalized_list);
ApplyCommandOnIntegrationRule(locate_hot_points, geom.getIntegRuleBasic(0), geom.begin("local"), geom.end("local"));
if(!delocalized_list.empty()) {
delocalized_list.sort();
delocalized_list.unique();
geom.updateDomains(delocalized_list);
delocalized_list.clear();
is_newly_delocalized=true;
......@@ -328,16 +326,8 @@ void TLSSolver::delocalizeMeanField() {
auto int_in_nonlocal=geom.getRegion("int_in_nonlocal");
lalg::xDenseMatrix coeffs(bnd_in_nonlocal.size(0), nonlocal.size(0));
// TODO refactor this it is too complicated
std::function<void (mVertex* v, double val)> set_val;
if(epsilon_ratio==0) {
set_val=[](mVertex* v, double val){};
}
else {
set_val=[this](mVertex* v, double val){ this->field.setVal(v, val); };
}
FastMarchingModeExtension([this](mVertex* v){ std::vector<double> vals; this->field.getVals(v, vals); return vals[0]; },
set_val,
[](mVertex* v, double val){},
nonlocal,
bnd_in_nonlocal.begin(0), bnd_in_nonlocal.end(0),
int_in_nonlocal.begin(0), int_in_nonlocal.end(0),
......
......@@ -31,6 +31,7 @@ void LocateHotPointsCommand::execute(xfem::xGeomElem* geom_integ) {
eval(geom_appro, geom_integ, val);
if (val>=threshold) {
list.push_back(geom_appro->getEntity());
break;
}
}
}
......
#!/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)