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

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

```
 1 2 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF 3 4 Primary Business: Queensland, Australia""" 5 __license__="""Licensed under the Open Software License version 3.0 6 7 from esys.escript import * 8 from esys.escript.linearPDEs import * 9 from esys import finley 10 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