Commit 10c05fbd authored by Alexis SALZMAN's avatar Alexis SALZMAN

Merge branch 'master' into para

Conflicts:
	CMakeLists.txt
parents 5b561468 a68a8e41
......@@ -104,7 +104,7 @@ if(FASTMARCHINGINTERFACE)
add_definitions(-DHAVE_FASTMARCHING)
find_package(FastMarching)
list(APPEND EXTERNAL_INCLUDES ${FASTMARCHING_INCLUDE_DIR})
#list(APPEND EXTERNAL_LIBRARIES ${FASTMARCHING_LIBRARIES})
# list(APPEND EXTERNAL_LIBRARIES ${FASTMARCHING_LIBRARIES})
set(CXXFLAGS "${CXXFLAGS} -std=c++11")
endif(FASTMARCHINGINTERFACE)
......
......@@ -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;
}
......
This diff is collapsed.
/*
This source code is subject to non-permissive licence,
see the TLSDuctile/LICENSE file for conditions.
*/
#ifndef _TLSDuctile_FormulationQSRemeshAniso_h_
#define _TLSDuctile_FormulationQSRemeshAniso_h_
#include "FormulationQS.h"
......@@ -28,6 +24,9 @@ public:
xfem::xLevelSet lsetAdapt;
std::function<void (xfem::xMesh&)> tagger;
xMesh *backupMesh;
private:
void updateLevelSetField();
void xReadFromParaview(std::string, std::string, xfem::xField&,xfem::xMesh&);
......@@ -38,8 +37,82 @@ private:
void projectPhiOnMesh(xMesh *mesh_new, std::map<mEntity *, double> &projectedValues);
void xAssignPhiValues(xField &field, xMesh *mesh, std::map<mEntity *, double> &projectedValues);
std::list<AOMD::mEntity*> constructDelocalizedList(xfem::xLevelSet&);
const double epsilon_ratio;
bool doRestart;
};
#include "mPoint.h"
template <class T>
class xEvalOnOtherMesh : public xfem::xEval<T>
{
public :
//Constructeur **********
xEvalOnOtherMesh(const xfem::xEval<T> &evalInit_, xMesh& meshInit_) : evalInit(evalInit_),meshInit(meshInit_) {}
//Operateur *****************
void operator()(const xGeomElem* geo_appro, const xGeomElem* geo_integ, T& resu) const
{
//On recupere les coordonnees globales
Trellis_Util::mPoint xyz = geo_integ->getXYZ(); //On recupere les coordonnees du point considere
//Recharge de l'element d'approx =============
std::set<AOMD::mEntity*> elts; //On cree un "set" d'elements
meshInit.locateElement(xyz,elts); //On cherche tous les elements qui contiennent xyz et on les met dans le set
//if(elts.begin()==elts.end()){std::cout<<"ERREUR dans le xEvalWiOnAll : Le noeud "<<xyz<<" n'est pas dans le maillage\nContenu du-dit maillage :\n";meshInit.printAll();throw;}
if(elts.begin()==elts.end()){
std::cout<<"ERREUR dans xEvalOnOtherMesh : Le noeud "<<xyz(0) <<" "<< xyz(1) <<" n'est pas dans le maillage"<<endl;
resu *= 0.;
return;
}
//on extrait l'element parmis tous les elements
AOMD::mEntity* eltInit=0; //On cree un element unique
eltInit=*(elts.begin()); //...et on dit que c'est le premier element du set
//Plus simple ici : les noeuds sont dupliques : xyz doit forcement etre place sur un noeud de meshInit
//donc pas la peine de chercher la partition (qui n'existe pas, vue la maniere dont fonctionne le code ici)
#if 0
//Recuperation de l'element d'integration (sous-element) =========
xfem::xPartition partition;
xMesh::getPartition(eltInit, partition, xfem::xAcceptAll());//Recupere les sous elements
mEntity* elemIntegInit=0; //Element selectionne
Trellis_Util::LagrangeMappingBuilder mapping_builder;//Mapping
double u,v,w;//??
for(xfem::xPartition::iterator it=partition.begin();it!=partition.end();it++)
{
mEntity* e=*it;
Trellis_Util::Mapping* mapping = mapping_builder.BuildMapping(e);
if (mapping->interiorCheck (e, xyz, u, v, w))
elemIntegInit=*it;
delete mapping;
}
xGeomElem geoApproInit(eltInit); //Espace d'approximation (l'element entier)
xGeomElem geoIntegInit(elemIntegInit);//Espace d'inegration (sous-element)
//MAJ des coordonnees du point de calcul sur geo_Fine !
geoApproInit.setUVWForXYZ(xyz);
geoIntegInit.setUVWForXYZ(xyz);
#else
xGeomElem geoApproInit(eltInit);
xGeomElem geoIntegInit(eltInit);
geoApproInit.setUVWForXYZ(xyz);
geoIntegInit.setUVWForXYZ(xyz);
#endif
evalInit(&geoApproInit,&geoIntegInit,resu); //On evalue
} //Fin de l'operateur
private :
const xfem::xEval<T> &evalInit;
xfem::xMesh &meshInit;
};
#endif
......@@ -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) {
......
......@@ -419,6 +419,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) {
......
......@@ -289,8 +289,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();
// PARA: must be communicate to all proc
......@@ -353,16 +351,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")