project conditional statement to print values makes code slow
6 months ago by
from dolfin import * import matplotlib.pyplot as plt E = 29e3 K = 500.0 sy = 36.0 Area = 0.1 print 'properties:', E,K,sy, Area mesh = IntervalMesh(60,0,60) #plot(mesh, interactive = True) V = FunctionSpace(mesh,'CG',1) P1 = FunctionSpace(mesh,'DG',0) u = TrialFunction(V) v = TestFunction(V) F = E*inner(grad(u),grad(v))*dx a = lhs(F) L = rhs(F) u = Function(V) class right(SubDomain): def inside(self,x,on_boundary): tol = 1e-10 return abs(x-60.0) < tol and on_boundary class left(SubDomain): def inside(self,x,on_boundary): tol = 1e-10 return abs(x-0.0) < tol and on_boundary Right = right() Left = left() uml = Expression('t', t = 0.0) bcl= DirichletBC(V, Constant(0.0), Left) bcr = DirichletBC(V, uml, Right) bc = [bcl, bcr] t = 0.0 dt = 0.005 num_steps = (0.5-0.0)/dt num_steps = 100 delta = 0.00005 pStrain = 0.0 alphan = 0.0 pStrainN = 0.0 def stre(Gamma,E,Tstress): return (1.0-(Gamma*E)/abs(Tstress))*Tstress def E_val(K,E): return E*K/(E+K) def alpha_val(alphan,Gamma): return alphan+Gamma def pStrain_val(pStrainN,Gamma,Tstress ): return pStrainN+ Gamma*sign(Tstress) for n in range(num_steps): t+=dt uml.t = t solve(a==L, u, bc) eps = u.dx(0) #print u.vector().norm('l2') ep = project(eps,P1) Tstrain = eps + delta Tstress = E*(Tstrain - pStrain) pStrain = pStrainN fun = abs(Tstress) -(sy + K*alphan) #fu = project(fun,P1) stress = Tstress st = project(E*(Tstrain),P1) right = Point(60.0) print 'str:',ep(right), 'Tstress', st(right) Gamma = fun/(E+K) stress = conditional(le(fun,0.0),Tstress,stre(Gamma,E,Tstress)) C_ep = conditional(le(fun,0.0),E, E_val(K,E) ) pStrain = conditional(le(fun,0.0),pStrain,pStrain_val(pStrainN,Gamma,Tstress ) ) alphaT = conditional(le(fun,0.0),alphan,alpha_val(alphan,Gamma)) alphan = alphaT pStraiN = pStrain if n>=20: stt = project(stress,P1) print 'stress:', stt(right)
Thanks for the help
Community: FEniCS Project
6 months ago by
Please login to add an answer/comment or follow this question.