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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 996 by gross, Tue Feb 20 06:19:30 2007 UTC revision 997 by gross, Mon Feb 26 06:31:45 2007 UTC
# Line 19  __version__="$Revision$" Line 19  __version__="$Revision$"
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):
# Line 49  class MiningHistory(Model): Line 51  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
Legend:
Removed from v.996  
changed lines
  Added in v.997
  ViewVC Help
Powered by ViewVC 1.1.26