# Contents of /trunk/finley/test/python/debugIntegrate.py

Revision 147 - (show annotations)
Fri Aug 12 01:45:47 2005 UTC (15 years ago) by jgs
Original Path: trunk/esys2/finley/test/python/debugIntegrate.py
File MIME type: text/x-python
File size: 1861 byte(s)
```erge of development branch dev-02 back to main trunk on 2005-08-12

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

## Properties

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