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

Annotation of /trunk/finley/test/python/debugIntegrate.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 617 - (hide annotations)
Wed Mar 22 02:58:17 2006 UTC (14 years, 4 months ago) by elspeth
File MIME type: text/x-python
File size: 2149 byte(s)
More copyright.

1 jgs 147
2 elspeth 617 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
3     http://www.access.edu.au
4     Primary Business: Queensland, Australia"""
5     __license__="""Licensed under the Open Software License version 3.0
6     http://www.opensource.org/licenses/osl-3.0.php"""
7 jgs 149 from esys.escript import *
8     from esys.escript.linearPDEs import *
9     from esys import finley
10 jgs 147
11     def makeParallelFaultSystem2D(xne=1,yne=1,xl=1.,yfaults=[0.,1.],integrationOrder=0,p1=0,p2=0):
12     """
13     generates a 2D rectangular domain and a system of parallel faults at yfaults
14    
15     input data :
16     xne - the number of elements over the horizontal edge of the rectangular mesh
17     xl - length of the horizontal edge of the domain
18     yfault[y0,...yi,...y1] where y0 and y1 are the extremity of the block and yi the fault position in the (Oy) direction
19     output data :
20     a mesh with horizontal discontinuities at yi positions
21     """
22    
23     #yne=xne*1./(xl+1.e-15)
24    
25     mshs=[]
26     for i in range(len(yfaults)-1):
27     yl=yfaults[i+1]-yfaults[i]
28     msh=finley.Rectangle(xne,yne,1,xl,yl,p1,p2,useElementsOnFace=True)
29     msh1=msh
30     print "mesh done"
31     msh1.setX(msh.getX()+[0.,yfaults[i]])
32     print "Xnew"
33     mshs.append(msh1)
34     print "apres append"
35     #results=finley.joinFaces(mshs)
36     mesh_joint=finley.JoinFaces(mshs)
37     print "apres faces jointes"
38     return mesh_joint
39    
40     def ContractProd(M1,M2):
41     msh=M1.getDomain()
42     ContractProd=inner(M1,M2)
43     return ContractProd
44    
45     def GetStrainEnergy(stressStateInit,uStateInit,stress,u):
46     strainStateInit=getStrain(uStateInit)
47     strain=getStrain(u)
48     E_Strain0=ContractProd(stressStateInit,strainStateInit).integrate()[0]
49     E_StrainCurrent=ContractProd(stress,strain).integrate()[0]
50     E_StrainE=-E_Strain0+E_StrainCurrent
51     return E_StrainE
52    
53     def main():
54     xne=10
55     yne=20
56     xl=10
57     yl=10
58     yfaults=[0.,5.,10.]
59     numberfaults=1
60    
61     mesh=makeParallelFaultSystem2D(xne,yne,xl,yfaults,0,0,0)
62    
63     t=Tensor(10.,Function(mesh))
64     t1=Tensor(20.,Function(mesh))
65    
66     contract_product=ContractProd(t,t1)
67     integral=integrate(contract_product)
68    
69     print "integral",integral
70    
71     s=Scalar(10.,Function(mesh))
72     print integrate(s)
73    
74     print "done"
75    
76     main()

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26