| 1 |
jgs |
147 |
|
| 2 |
elspeth |
617 |
__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 |
jgs |
149 |
from esys.escript import * |
| 8 |
|
|
from esys.escript.linearPDEs import * |
| 9 |
|
|
from esys import finley |
| 10 |
jgs |
147 |
|
| 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() |
Parent Directory
| 
Revision Log