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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1809 - (show annotations)
Thu Sep 25 06:43:44 2008 UTC (10 years, 11 months ago) by ksteube
File MIME type: text/x-python
File size: 2625 byte(s)
Copyright updated in all python files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26