1 |
|
2 |
__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 |
from esys.escript import * |
8 |
from esys.escript.linearPDEs import * |
9 |
from esysi 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 main(): |
41 |
|
42 |
xne=10 |
43 |
yne=20 |
44 |
xl=10 |
45 |
yl=10 |
46 |
yfaults=[0.,5.,10.] |
47 |
numberfaults=1 |
48 |
|
49 |
mesh=makeParallelFaultSystem2D(xne,yne,xl,yfaults,0,0,0) |
50 |
|
51 |
v=Vector(10.,ContinuousFunction(mesh)) |
52 |
|
53 |
on_elem_v=interpolate(v,Function(mesh)) |
54 |
|
55 |
value_integratedEL_vec=integrate(on_elem_v) |
56 |
|
57 |
print "vector--------" |
58 |
print "value interpolation",value_integratedEL_vec |
59 |
print "--------" |
60 |
|
61 |
main() |