/[escript]/trunk/finley/test/python/linearElastic.py
ViewVC logotype

Diff of /trunk/finley/test/python/linearElastic.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3410 by caltinay, Fri Nov 12 01:19:02 2010 UTC revision 3411 by gross, Thu Dec 9 22:42:02 2010 UTC
# Line 24  from esys.escript.linearPDEs import Line Line 24  from esys.escript.linearPDEs import Line
24  from esys import finley  from esys import finley
25  from esys.weipa import saveVTK  from esys.weipa import saveVTK
26    
27  press0=1.  pres0=-100.
28  lamb=1.  lame=1.
29  nu=0.3  mu=0.3
30    rho=1.
31    g=9.81
32    
33  # this sets the hook tensor:  # generate mesh: here 20x20 mesh of order 1
34  def setHookTensor(w,l,n):  domain=finley.Rectangle(20,20,1,l0=1.0,l1=1.0)
    C=Tensor4(0.,w)  
    for i in range(w.getDim()):  
      for j in range(w.getDim()):  
           C[i,i,j,j]+=l  
           C[j,i,j,i]+=n  
           C[j,i,i,j]+=n  
   
    return C  
     
 # generate mesh: here 10x20 mesh of order 2  
 domain=finley.Rectangle(10,20,1,l0=0.5,l1=1.0)  
 # get handel to nodes and elements:  
 e=Function(domain)  
 fe=FunctionOnBoundary(domain)  
 n=ContinuousFunction(domain)  
35  #  #
36  # set a mask msk of type vector which is one for nodes and components set be a constraint:  # set a mask msk of type vector which is one for nodes and components set be a constraint:
37  #  #
38  msk=whereZero(n.getX()[0])*[1.,1.]  msk=whereZero(domain.getX()[0])*[1.,1.]
39  #  #
40  #  set the normal stress components on face elements.  #  set the normal stress components on face elements.
41  #  faces tagged with 21 get the normal stress [0,-press0].  #  faces tagged with 21 get the normal stress [0,-press0].
42  #  #
43  # now the pressure is set to zero for x0 coordinates bigger then 0.1  # now the pressure is set to zero for x0 coordinates equal 1. (= right face)
44  press=whereNegative(fe.getX()[0]-0.1)*200000.*[1.,0.]  press=whereZero(FunctionOnBoundary(domain).getX()[0]-1.)*pres0*[1.,0.]
45  # assemble the linear system:  # assemble the linear system:
46  mypde=LinearPDE(domain)  mypde=LinearPDE(domain)
47  mypde.setValue(A=setHookTensor(e,lamb,nu),y=press,q=msk,r=[0,0])  k3=kronecker(domain)
48    k3Xk3=outer(k3,k3)
49    
50    mypde.setValue(A=mu * ( swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3) ) + lame*k3Xk3,
51                   Y=[0,-g*rho],
52                   y=press,
53                   q=msk,r=[0,0])
54  mypde.setSymmetryOn()  mypde.setSymmetryOn()
55  mypde.getSolverOptions().setVerbosityOn()  mypde.getSolverOptions().setVerbosityOn()
56  mypde.getSolverOptions().setPreconditioner(mypde.getSolverOptions().AMG)  # use direct solver (default is iterative)
57    mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().DIRECT)
58    # mypde.getSolverOptions().setPreconditioner(mypde.getSolverOptions().AMG)
59  # solve for the displacements:  # solve for the displacements:
60  u_d=mypde.getSolution()  u_d=mypde.getSolution()
61  # get the gradient and calculate the stress:  # get the gradient and calculate the stress:
62  g=grad(u_d)  g=grad(u_d)
63  stress=lamb*trace(g)*kronecker(domain)+nu*(g+transpose(g))  stress=lame*trace(g)*kronecker(domain)+mu*(g+transpose(g))
64  # write the hydrostatic pressure:  # write the hydrostatic pressure:
65  saveVTK("result.vtu",displacement=u_d,pressure=trace(stress)/domain.getDim())  saveVTK("result.vtu",displacement=u_d,pressure=trace(stress)/domain.getDim())

Legend:
Removed from v.3410  
changed lines
  Added in v.3411

  ViewVC Help
Powered by ViewVC 1.1.26