| 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 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() |
Parent Directory
| 
Revision Log