19 |
__date__="$Date$" |
__date__="$Date$" |
20 |
|
|
21 |
from esys.escript.modelframe import Model |
from esys.escript.modelframe import Model |
22 |
|
from esys.escript import * |
23 |
|
from esys.escript.linearPDEs import LinearPDE |
24 |
from mines import parse |
from mines import parse |
25 |
|
|
26 |
class MiningHistory(Model): |
class MiningHistory(Model): |
51 |
|
|
52 |
def doStepPreprocessing(self, dt): |
def doStepPreprocessing(self, dt): |
53 |
self.mass_changes=self.__minesite.getMassChanges(self.t) |
self.mass_changes=self.__minesite.getMassChanges(self.t) |
|
print self.t,":",self.mass_changes |
|
54 |
|
|
55 |
def getSafeTimeStepSize(self, dt): |
def getSafeTimeStepSize(self, dt): |
|
print "new time marker:",self.t,":",self.__minesite.getNextTimeMarker(self.t) |
|
56 |
return self.__minesite.getNextTimeMarker(self.t)-self.t |
return self.__minesite.getNextTimeMarker(self.t)-self.t |
57 |
|
|
58 |
|
|
59 |
|
class DensityChange(Model): |
60 |
|
""" |
61 |
|
translates local mass extraction into density changes. |
62 |
|
"local" refers a taged region. |
63 |
|
|
64 |
|
@ivar domain: mining region |
65 |
|
@type domain: L{Domian} |
66 |
|
@ivar tag_map: a tagmap |
67 |
|
@type tag_map: L{TagMap} |
68 |
|
@ivar mass_rate: rate of change mass in a tagged region |
69 |
|
@type mass_rate: C{dict} |
70 |
|
@ivar density_rate: density in each region |
71 |
|
@type density_rate: C{Data} |
72 |
|
""" |
73 |
|
def __init__(self,**kwargs): |
74 |
|
""" |
75 |
|
""" |
76 |
|
super(DensityChange,self).__init__(**kwargs) |
77 |
|
self.declareParameter(domain=None, |
78 |
|
tag_map=None, |
79 |
|
mass_rate={}, |
80 |
|
density_rate=None) |
81 |
|
|
82 |
|
def doInitialization(self): |
83 |
|
""" |
84 |
|
initialize time integration |
85 |
|
""" |
86 |
|
self.__volumes={} |
87 |
|
for i in self.tag_map.getNames(): |
88 |
|
c=Scalar(0.,Function(self.domain)) |
89 |
|
self.tag_map.insert(c,**{i : 1.}) |
90 |
|
c=integrate(c) |
91 |
|
if c>0: self.__volumes[i] = c |
92 |
|
|
93 |
|
def doStepPreprocessing(self, dt): |
94 |
|
self.density_rate=Scalar(0.,Function(self.domain)) |
95 |
|
d={} |
96 |
|
for i in self.__volumes: |
97 |
|
if self.mass_rate.has_key(i): d[i]=self.mass_rate[i]/self.__volumes[i] |
98 |
|
self.tag_map.insert(self.density_rate,**d) |
99 |
|
|
100 |
|
class LinearElasticStressChange(Model): |
101 |
|
""" |
102 |
|
calculates the stress according to an initial garvity field and a consecutive |
103 |
|
chenge of density based on linar elastic material model. It is assumed that |
104 |
|
the lame coefficients don't change over time, the direction of gravity is |
105 |
|
pointing into the x_{dim} direction. |
106 |
|
|
107 |
|
@note: add stress changes due to tectonic loading and slip |
108 |
|
|
109 |
|
@ivar domain: mining region |
110 |
|
@type domain: L{Domian} |
111 |
|
@ivar tag_map: a tagmap |
112 |
|
@type tag_map: L{TagMap} |
113 |
|
@ivar displacement: displacement field |
114 |
|
@type displacement: L{Vector} or C{None} |
115 |
|
@ivar stress: displacement field |
116 |
|
@type stress: L{Vector} or C{None} |
117 |
|
@ivar density: initial density distribution |
118 |
|
@type density: C{float} or L{Scalar} |
119 |
|
@ivar density_rate: density rate by tag (may be changed of time) |
120 |
|
@type density_rate: C{dict} |
121 |
|
@ivar lame_lambda: elasticity coefficient lambda (assumed to be constant over time) |
122 |
|
@type lame_lambda: C{float} or L{Scalar} |
123 |
|
@ivar lame_mu: elasticity coefficient mu (assumed to be constant over time) |
124 |
|
@type lame_mu: C{float} or L{Scalar} |
125 |
|
@ivar location_of_fixed_displacement: mask of locations and component with zero displacements |
126 |
|
@type location_of_fixed_displacement: L{Vector} or C{None} |
127 |
|
""" |
128 |
|
def __init__(self,**kwargs): |
129 |
|
""" |
130 |
|
""" |
131 |
|
super(LinearElasticStressChange,self).__init__(**kwargs) |
132 |
|
self.declareParameter(domain=None, |
133 |
|
tag_map=None, |
134 |
|
displacement=None, \ |
135 |
|
stress=None, \ |
136 |
|
density=1., \ |
137 |
|
lame_lambda=2., \ |
138 |
|
lame_mu=1., \ |
139 |
|
location_of_fixed_displacement=None, \ |
140 |
|
density_rate=None) |
141 |
|
|
142 |
|
def doInitialization(self): |
143 |
|
""" |
144 |
|
initialize time integration |
145 |
|
""" |
146 |
|
if self.stress == None: self.stress=Tensor(0,Function(self.domain)) |
147 |
|
if self.displacement==None: self.displacement=Scalar(0,Solution(self.domain)) |
148 |
|
self.__pde=LinearPDE(self.domain) |
149 |
|
self.__pde.setSymmetryOn() |
150 |
|
A =Tensor4(0.,Function(self.domain)) |
151 |
|
for i in range(self.domain.getDim()): |
152 |
|
for j in range(self.domain.getDim()): |
153 |
|
A[i,j,j,i] += self.lame_mu |
154 |
|
A[i,j,i,j] += self.lame_mu |
155 |
|
A[i,i,j,j] += self.lame_lambda |
156 |
|
self.__pde.setValue(A=A,q=self.location_of_fixed_displacement) |
157 |
|
|
158 |
|
self.__stress=self.stress |
159 |
|
self.__displacement=self.displacement |
160 |
|
self.__first=True |
161 |
|
|
162 |
|
def doInitialStep(self): |
163 |
|
""" |
164 |
|
""" |
165 |
|
self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]*9.81*self.density) |
166 |
|
ddisp=self.__pde.getSolution() |
167 |
|
deps=symmetric(grad(ddisp)) |
168 |
|
self.stress=self.stress+self.lame_mu*deps+self.lame_lambda*trace(deps)*kronecker(Function(self.domain)) |
169 |
|
self.displacement=self.displacement+ddisp |
170 |
|
self.__first=False |
171 |
|
|
172 |
|
def terminateInitialIteration(self): |
173 |
|
return not self.__first |
174 |
|
|
175 |
|
def terminateIteration(self): |
176 |
|
return not self.__first |
177 |
|
|
178 |
|
def doInitialStepPostprocessing(self): |
179 |
|
self.__stress=self.stress |
180 |
|
self.__displacement=self.displacement |
181 |
|
|
182 |
|
def doStepPreprocessing(self,dt): |
183 |
|
self.stress=self.__stress |
184 |
|
self.displacement=self.__displacement |
185 |
|
self.__first=True |
186 |
|
|
187 |
|
def doStep(self): |
188 |
|
""" |
189 |
|
""" |
190 |
|
if not self.density_rate == None: |
191 |
|
self.__pde.setValue(Y=-kronecker(Function(self.domain))[d-1]*9.81*self.density_rate) |
192 |
|
ddisp=self.__pde.getSolution() |
193 |
|
deps=symmetric(grad(ddisp)) |
194 |
|
self.stress=self.stress+dt*(self.lame_mu*deps+self.lame_lambda*trace(deps)*kronecker(Function(self.domain))) |
195 |
|
self.displacement=self.displacement+dt*ddisp |
196 |
|
self.__first=False |
197 |
|
|
198 |
|
def doStepPostprocessing(self, dt): |
199 |
|
self.__stress=self.stress |
200 |
|
self.__displacement=self.displacement |