test/python/debugIntegrate.py


from escript.escript import *
from escript.linearPDEs import *
from finley import finley

def makeParallelFaultSystem2D(xne=1,yne=1,xl=1.,yfaults=[0.,1.],integrationOrder=0,p1=0,p2=0):
 """
 generates a 2D rectangular domain and a system of parallel faults at yfaults 
   
 input data :
 xne - the number of elements over the horizontal edge of the rectangular mesh
 xl  - length of the horizontal edge of the domain
 yfault[y0,...yi,...y1] where y0 and y1 are the extremity of the block and yi the fault position in the (Oy) direction
   output data :
        a mesh with horizontal discontinuities at yi positions
 """

 #yne=xne*1./(xl+1.e-15)
 
 mshs=[]
 for i in range(len(yfaults)-1):
     yl=yfaults[i+1]-yfaults[i]
     msh=finley.Rectangle(xne,yne,1,xl,yl,p1,p2,useElementsOnFace=True)
     msh1=msh
     print "mesh done"
     msh1.setX(msh.getX()+[0.,yfaults[i]])
     print "Xnew"
     mshs.append(msh1)
     print "apres append"
     #results=finley.joinFaces(mshs)
     mesh_joint=finley.JoinFaces(mshs)
     print "apres faces jointes"
 return mesh_joint

def ContractProd(M1,M2):
 msh=M1.getDomain()
 ContractProd=inner(M1,M2)
 return ContractProd

def GetStrainEnergy(stressStateInit,uStateInit,stress,u):
 strainStateInit=getStrain(uStateInit)
 strain=getStrain(u)
 E_Strain0=ContractProd(stressStateInit,strainStateInit).integrate()[0]
 E_StrainCurrent=ContractProd(stress,strain).integrate()[0]
 E_StrainE=-E_Strain0+E_StrainCurrent
 return E_StrainE

def main():
 xne=10
 yne=20
 xl=10
 yl=10
 yfaults=[0.,5.,10.]
 numberfaults=1

 mesh=makeParallelFaultSystem2D(xne,yne,xl,yfaults,0,0,0)

 t=Tensor(10.,Function(mesh))
 t1=Tensor(20.,Function(mesh))

 contract_product=ContractProd(t,t1)
 integral=integrate(contract_product)
   
 print "integral",integral

 s=Scalar(10.,Function(mesh))
 print integrate(s)
   
 print  "done"
   
main()
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()
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision