/[escript]/trunk/modellib/py_src/crustal/setups.py
ViewVC logotype

Contents of /trunk/modellib/py_src/crustal/setups.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1044 - (show annotations)
Mon Mar 19 07:29:31 2007 UTC (12 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 9155 byte(s)
clear name tagging is supported now.
1 """
2 set up of mining activities in modelframe
3
4 @var __author__: name of author
5 @var __licence__: licence agreement
6 @var __url__: url entry point on documentation
7 @var __version__: version
8 @var __date__: date of the version
9 """
10
11 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
12 http://www.access.edu.au
13 Primary Business: Queensland, Australia"""
14 __license__="""Licensed under the Open Software License version 3.0
15 http://www.opensource.org/licenses/osl-3.0.php"""
16 __author__="Lutz Gross, l.gross@uq.edu.au"
17 __url__="http://www.iservo.edu.au/esys/escript"
18 __version__="$Revision$"
19 __date__="$Date$"
20
21 from esys.escript.modelframe import Model, ParameterSet
22 from esys.escript import *
23 from esys.escript.linearPDEs import LinearPDE
24 from mines import parse
25 import numarray
26
27 class MiningHistory(Model):
28 """
29 manages the history of mines. It mandels the time steping according to the available
30 data and a dictionary of mass changes per year for all the mines.
31
32 @ivar history: mine site history file.
33 @type history: C{DataSource}
34 @ivar mass_changes: current mass change per year.
35 @type mass_changes: C{DataSource}
36
37 """
38 def __init__(self,**kwargs):
39 """
40 """
41 super(MiningHistory,self).__init__(**kwargs)
42 self.declareParameter(t=0.,
43 history=None,
44 mass_changes=None)
45
46 def doInitialization(self):
47 """
48 initialize time integration
49 """
50 self.__minesite=parse(open(self.history.getLocalFileName(),'r'))
51 self.mass_changes=self.__minesite.getMassChanges(self.t)
52
53 def doStepPreprocessing(self, dt):
54 self.mass_changes=self.__minesite.getMassChanges(self.t)
55
56 def getSafeTimeStepSize(self, dt):
57 return self.__minesite.getNextTimeMarker(self.t)-self.t
58
59
60 class DensityChange(Model):
61 """
62 translates local mass extraction into density changes.
63 "local" refers a taged region.
64
65 @ivar domain: mining region
66 @type domain: L{Domian}
67 @ivar mass_rate: rate of change mass in a tagged region
68 @type mass_rate: C{dict}
69 @ivar density_rate: density in each region
70 @type density_rate: C{Data}
71 """
72 def __init__(self,**kwargs):
73 """
74 """
75 super(DensityChange,self).__init__(**kwargs)
76 self.declareParameter(domain=None,
77 mass_rate={},
78 density_rate=None)
79
80 def doInitialization(self):
81 """
82 initialize time integration
83 """
84 self.__volumes={}
85 for i in getTagNames(self.domain):
86 c=integrate(insertTaggedValues(Function(self.domain),**{i : 1.}))
87 if c>0: self.__volumes[i] = c
88
89 def doStepPreprocessing(self, dt):
90 d={}
91 for i in self.__volumes:
92 if self.mass_rate.has_key(i): d[i]=-self.mass_rate[i]/self.__volumes[i]
93 self.density_rate=inserTaggedValues(Scalar(0.,Function(self.domain)),**d)
94
95 class LinearElasticStressChange(Model):
96 """
97 calculates the stress according to an initial garvity field and a consecutive
98 chenge of density based on linar elastic material model. It is assumed that
99 the lame coefficients don't change over time, the direction of gravity is
100 pointing into the x_{dim} direction.
101
102 @note: add stress changes due to tectonic loading and slip
103
104 @ivar domain: mining region
105 @type domain: L{Domian}
106 @ivar displacement: displacement field
107 @type displacement: L{Vector} or C{None}
108 @ivar stress: displacement field
109 @type stress: L{Vector} or C{None}
110 @ivar density: initial density distribution
111 @type density: C{float} or L{Scalar}
112 @ivar density_rate: density rate by tag (may be changed of time)
113 @type density_rate: C{dict}
114 @ivar lame_lambda: elasticity coefficient lambda (assumed to be constant over time)
115 @type lame_lambda: C{float} or L{Scalar}
116 @ivar lame_mu: elasticity coefficient mu (assumed to be constant over time)
117 @type lame_mu: C{float} or L{Scalar}
118 @ivar location_of_fixed_displacement: mask of locations and component with zero displacements
119 @type location_of_fixed_displacement: L{Vector} or C{None}
120 """
121 def __init__(self,**kwargs):
122 """
123 """
124 super(LinearElasticStressChange,self).__init__(**kwargs)
125 self.declareParameter(domain=None,
126 displacement=None, \
127 stress=None, \
128 density=1., \
129 lame_lambda=2., \
130 lame_mu=1., \
131 location_of_fixed_displacement=None, \
132 density_rate=None)
133
134 def doInitialization(self):
135 """
136 initialize time integration
137 """
138 if self.stress == None: self.stress=Tensor(0,Function(self.domain))
139 if self.displacement==None: self.displacement=Scalar(0,Solution(self.domain))
140 self.__pde=LinearPDE(self.domain)
141 self.__pde.setSymmetryOn()
142 A =Tensor4(0.,Function(self.domain))
143 for i in range(self.domain.getDim()):
144 for j in range(self.domain.getDim()):
145 A[i,j,j,i] += self.lame_mu
146 A[i,j,i,j] += self.lame_mu
147 A[i,i,j,j] += self.lame_lambda
148 self.__pde.setValue(A=A,q=self.location_of_fixed_displacement)
149
150 self.__stress=self.stress
151 self.__displacement=self.displacement
152 self.__first=True
153
154 def doInitialStep(self):
155 """
156 """
157 d=self.domain.getDim()
158 self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]*9.81*self.density)
159 ddisp=self.__pde.getSolution()
160 deps=symmetric(grad(ddisp))
161 self.stress=self.stress+self.lame_mu*deps+self.lame_lambda*trace(deps)*kronecker(Function(self.domain))
162 self.displacement=self.displacement+ddisp
163 self.__first=False
164
165 def terminateInitialIteration(self):
166 return not self.__first
167
168 def terminateIteration(self):
169 return not self.__first
170
171 def doInitialPostprocessing(self):
172 self.__stress=self.stress
173 self.__displacement=self.displacement
174 self.trace("initial displacement range : %s: %s"%(inf(self.displacement[2]),sup(self.displacement[2])))
175 self.trace("initial stress range : %s: %s"%(inf(self.stress),sup(self.stress)))
176
177 def doStepPreprocessing(self,dt):
178 self.stress=self.__stress
179 self.displacement=self.__displacement
180 self.__first=True
181
182 def doStep(self,dt):
183 """
184 """
185 d=self.domain.getDim()
186 if not self.density_rate == None:
187 self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]*9.81*self.density_rate)
188 ddisp=self.__pde.getSolution()
189 deps=symmetric(grad(ddisp))
190 self.stress=self.stress+dt*(self.lame_mu*deps+self.lame_lambda*trace(deps)*kronecker(Function(self.domain)))
191 self.displacement=self.displacement+dt*ddisp
192 self.__first=False
193
194 def doStepPostprocessing(self, dt):
195 self.__stress=self.stress
196 self.__displacement=self.displacement
197 self.trace("displacement range : %s: %s"%(inf(self.displacement[2]),sup(self.displacement[2])))
198 self.trace("stress range : %s: %s"%(inf(self.stress),sup(self.stress)))
199
200 class CoulombFailureStress(ParameterSet):
201 """
202 calculates the Coulomb failure stress an planes of a given orientation.
203
204 @ivar domain: mining region
205 @type domain: L{Domian}
206 @ivar displacement: displacement field
207 @type displacement: L{Vector} or C{None}
208 @ivar stress: displacement field
209 @type stress: L{Vector} or C{None}
210 @ivar density: initial density distribution
211 @type density: C{float} or L{Scalar}
212 @ivar density_rate: density rate by tag (may be changed of time)
213 @type density_rate: C{dict}
214 @ivar lame_lambda: elasticity coefficient lambda (assumed to be constant over time)
215 @type lame_lambda: C{float} or L{Scalar}
216 @ivar lame_mu: elasticity coefficient mu (assumed to be constant over time)
217 @type lame_mu: C{float} or L{Scalar}
218 @ivar location_of_fixed_displacement: mask of locations and component with zero displacements
219 @type location_of_fixed_displacement: L{Vector} or C{None}
220 """
221 def __init__(self,**kwargs):
222 """
223 """
224 super(CoulombFailureStress,self).__init__(**kwargs)
225 self.declareParameter(stress=numarray.zeros((3,3)),
226 friction_coefficient=0.,
227 normal=numarray.array([1.,0.,0.]))
228 def cfs(self):
229 """
230 returns Coulomb failure stress
231 """
232 sn=matrixmult(self.stress,self.normal)
233 nsn=inner(self.normal,sn)
234 nssn=inner(sn,sn)
235 return (sqrt(nssn-nsn**2)-self.friction_coefficient*nsn)/length(self.normal)
236

  ViewVC Help
Powered by ViewVC 1.1.26