/[escript]/trunk/escript/py_src/faultsystems.py
ViewVC logotype

Contents of /trunk/escript/py_src/faultsystems.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2647 - (show annotations)
Fri Sep 4 05:25:25 2009 UTC (9 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 14772 byte(s)
fault system add. There is still an example for the usage missing.
1 ########################################################
2 #
3 # Copyright (c) 2003-2009 by University of Queensland
4 # Earth Systems Science Computational Center (ESSCC)
5 # http://www.uq.edu.au/esscc
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 ########################################################
12
13 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
14 Earth Systems Science Computational Center (ESSCC)
15 http://www.uq.edu.au/esscc
16 Primary Business: Queensland, Australia"""
17 __license__="""Licensed under the Open Software License version 3.0
18 http://www.opensource.org/licenses/osl-3.0.php"""
19 __url__="https://launchpad.net/escript-finley"
20
21 from esys.escript import *
22 import numpy
23 import math
24
25 class FaultSystem:
26 """
27 The FaultSystem class defines a system of faults in the Earth crust.
28
29 A fault system is defined by set of faults index by a tag. Each fault is defined as a polygon describing the top of
30 the fault and - in case of a 3D model - a polygon describing the bottom of the fault. This class provides a mechanism
31 to parametrise each fault with the domain [0,w0_max] x [0, w1_max] (to [0,w0_max] in the 2D case).
32 """
33 NOTAG="__NOTAG__"
34 def __init__(self,dim=3):
35 """
36 Sets up the fault system
37
38 :param dim: spatial dimension
39 :type dim: ``int`` of value 2 or 3
40 """
41 if not (dim == 2 or dim == 3):
42 raise ValueError,"only dimension2 2 and 3 are supported."
43 self.__dim=dim
44 self.__faults={}
45 self.__length={}
46 self.__depth={}
47 self.__offsets={}
48 self.__w1_max={}
49 self.__w0_max={}
50 self.__center=None
51 self.__orientation = None
52
53 def getDim(self):
54 """
55 returns the spatial dimension
56 :rtype: ``int``
57 """
58 return self.__dim
59 def getLength(self,tag=None):
60 """
61 returns the unrolled length of fault ``tag``
62 :rtype: ``float``
63 """
64 if tag==None: tag=self.NOTAG
65 return self.__length[tag]
66
67 def getDepth(self,tag=None):
68 """
69 returns the medium depth of fault ``tag``
70 :rtype: ``float``
71 """
72 if tag==None: tag=self.NOTAG
73 return self.__depth[tag]
74
75 def getTags(self):
76 """
77 returns a list of the tags used by the fault system
78 :rtype: ``list``
79 """
80 return self.__faults.keys()
81
82 def getW0Range(self,tag=None):
83 """
84 returns the range of the parameterization in ``w0``
85 :rtype: two ``float``
86 """
87 return self.getW0Offsets(tag)[0], self.getW0Offsets(tag)[-1]
88
89 def getW1Range(self,tag=None):
90 """
91 returns the range of the parameterization in ``w1``
92 :rtype: two ``float``
93 """
94 if tag==None: tag=self.NOTAG
95 return -self.__w1_max[tag],0
96
97 def getW0Offsets(self, tag=None):
98 """
99 returns the offsets for the parametrization of fault ``tag``.
100
101 :return: the offsets in the parametrization
102 :rtype: ``list`` of ``float``
103 """
104 if tag==None: tag=self.NOTAG
105 return self.__offsets[tag]
106
107
108 def getFaultSegments(self, tag=None):
109 """
110 returns the polygons used to describe fault tagged by ``tag``
111
112 :param tag: the tag of the fault
113 :type tag: ``float`` or ``str``
114 :return: the list of vertices defining the top and the list of the vertices defining the bottom of the fault indexed by ``tag``. The coordinates are `numpy.ndarray'.
115 """
116 if tag==None: tag=self.NOTAG
117 return self.__faults[tag]
118
119 def getCenterOnSurface(self):
120 """
121 returns the center point of the fault system at the surface
122 :rtype: `numpy.ndarray`
123 """
124 if self.__center == None:
125 self.__center=numpy.zeros((3,),numpy.float)
126 counter=0
127 for t in self.getTags():
128 for s in self.getFaultSegments(t)[0]:
129 self.__center[:2]+=s[:2]
130 counter+=1
131 self.__center/=counter
132 return self.__center[:self.getDim()]
133
134 def getOrientationOnSurface(self):
135 """
136 returns the orientation of the fault system in RAD on the surface around the fault system center
137 :rtype: ``float``
138 """
139 if self.__orientation == None:
140 center=self.getCenterOnSurface()
141 covariant=numpy.zeros((2,2))
142 for t in self.getTags():
143 for s in self.getFaultSegments(t)[0]:
144 covariant[0,0]+=(center[0]-s[0])**2
145 covariant[0,1]+=(center[1]-s[1])*(center[0]-s[0])
146 covariant[1,1]+=(center[1]-s[1])**2
147 covariant[1,0]+=(center[1]-s[1])*(center[0]-s[0])
148 e, V=numpy.linalg.eigh(covariant)
149 if e[0]>e[1]:
150 d=V[:,1]
151 else:
152 d=V[:,0]
153 if abs(d[0])>0.:
154 self.__orientation=atan(d[1]/d[0])
155 else:
156 self.__orientation=math.pi/2
157 return self.__orientation
158 def transform(self, rot=0, shift=numpy.zeros((3,))):
159 """
160 applies a shift and a consecutive rotation in the x0x1 plane.
161
162 :param rot: rotation angle in RAD
163 :type rot: ``float``
164 :param shift: shift vector to be applied before rotation
165 :type shift: `numpy.ndarray` of size 2 or 3
166 """
167 if self.getDim() == 2:
168 mat=numpy.array([[cos(rot), -sin(rot)], [sin(rot), cos(rot)] ])
169 else:
170 mat=numpy.array([[cos(rot), -sin(rot),0.], [sin(rot), cos(rot),0.], [0.,0.,1.] ])
171 for t in self.getTags():
172 top=[]
173 for s in self.getFaultSegments(t)[0]: top.append(numpy.dot(mat,(s+shift[:self.getDim()])))
174 if self.getDim() == 3:
175 bottom=[]
176 for s in self.getFaultSegments(t)[1]: bottom.append(numpy.dot(mat,(s+shift[:self.getDim()])))
177 else:
178 bottom=None
179 self.addFault(top=top, tag=t, bottom=bottom, w0_offsets=self.getW0Offsets(t), w1_max=-self.getW1Range(t)[0])
180
181 def addFault(self, top, tag=None, bottom=None, w0_offsets=None, w1_max=None):
182 """
183 adds a new fault to the fault system. The fault is named by ``tag``.
184
185 a fault is subdivided into segments where segment ``i`` is represented by the four corner points
186 ``v0``=``top[i]``, ``v1``=``bottom[i]``, ``v2``=``bottom[i+1]`` and ``v3``==``top[i+1]``. The segment is parametrized
187 by ``w0`` and ``w1`` with ``w0_offsets[i]<=w0<=w0_offsets[i+1]`` and ``-w1_max<=w1<=0``. Moreover
188
189 - ``(w0,w1)=(w0_offsets[i] , 0)->v0``
190 - ``(w0,w1)=(w0_offsets[i] , -w1_max)->v1``
191 - ``(w0,w1)=(w0_offsets[i+1], -w1_max)->v2``
192 - ``(w0,w1)=(w0_offsets[i+1], 0)->v3``
193
194 If no ``w0_offsets`` is given,
195
196 - ``w0_offsets[0]=0``
197 - ``w0_offsets[i]=w0_offsets[i-1]+length(v0-v4)``
198
199 If no ``w1_max`` is given, the average fault depth (=length(v0-v1)) is used.
200
201 :param top: list of coordinates. top[i] and top[i+1] are the top corners of the i-th segment of the fault
202 :type top: ``list`` of objects which can be converted to `numpy.ndarray` vector of length 3 (or 2 in the 2D case)
203 :param tag: the tag of the fault. If fault ``tag`` already exists it is overwritten.
204 :type tag: ``float`` or ``str``
205 :param bottom: list of coordinates. bottom[i] and bottom[i+1] are the bottom corners of the i-th segment of the fault
206 :type bottom: ``list`` of objects which can be converted to `numpy.ndarray` vector of length 3 (or 2 in the 2D case)
207 :param w0_offsets: ``w0_offsets[i]`` defines the offset of the segment ``i`` in the fault to be used in the parametrization of the fault. If not present the cumulative length of the fault segments is used.
208 :type w0_offsets: ``list`` of ``float`` or ``None``
209 :param w1_max: the maximum value used for parametrization of the fault in the depth direction. If not present the mean depth of the fault segments is used.
210 :type w1_max: ``float``
211 :note: In the three dimensional case the lists ``bottom`` and ``top`` must have the same length.
212 """
213 if tag==None:
214 tag=self.NOTAG
215 else:
216 if self.NOTAG in self.getTags():
217 raise ValueError,'Attempt to add a fault with no tag to a set of existing faults'
218 if len(top)<2:
219 raise Value,"at least two vertices must be given."
220 if self.getDim()==2 and not bottom==None:
221 raise ValueError,'Spatial dimension two does not support bottom polygon for faults.'
222 if not bottom == None:
223 if len(top) != len(bottom):
224 raise ValueError,'length of top and bottom polygon must be identical.'
225 if w0_offsets != None:
226 if len(w0_offsets) != len(top):
227 raise ValueError,'expected length of w0_offsets is %s'%(len(top))
228 self.__center=None
229 self.__orientation = None
230 #
231 # translate to numpy
232 #
233 if self.getDim()==2:
234 self.__faults[tag]=([ numpy.array(t[:2],numpy.double) for t in top ],)
235 else:
236 self.__faults[tag]=( [ numpy.array(t[:3],numpy.double) for t in top ], [ numpy.array(b[:3],numpy.double) for b in bottom ])
237 faults=self.__faults[tag]
238 len_faults=len(faults[0])
239 #
240 # calculate the depth
241 #
242 if self.getDim()==2:
243 self.__depth[tag]=0.
244 else:
245 self.__depth[tag]=sum([numpy.linalg.norm(faults[1][i]-faults[0][i]) for i in xrange(len_faults)])/len_faults
246 if w1_max==None or self.getDim()==2: w1_max=self.__depth[tag]
247 self.__length[tag]=sum([numpy.linalg.norm(faults[0][i]-faults[0][i-1]) for i in xrange(1,len_faults)])
248 self.__w1_max[tag]=w1_max
249 #
250 #
251 #
252 if w0_offsets!=None:
253 self.__offsets[tag]=w0_offsets
254 else:
255 self.__offsets[tag]=[0.]
256 for i in xrange(1,len_faults):
257 self.__offsets[tag].append(self.__offsets[tag][-1]+numpy.linalg.norm(faults[0][i]-faults[0][i-1]))
258 w0_max=max(self.__offsets[tag])
259 self.__w0_max[tag]=w0_max
260
261 def getMaxValue(self,f, tol=sqrt(EPSILON)):
262 """
263 returns the maximum value of ``f``, the fault and the location on the fault in fault coordinates.
264
265 :param f: a distribution of values
266 :type f: `escript.Data`
267 :param tol: relative tolerance used to decide if point is on fault
268 :type f: ``tol``
269 :return: the maximum value of across all faults in the fault system, the fault tag the maximum is taken, and the coordinates of the location in the coordinates of the fault. The returned fault tag is ``None`` if no points on the fault are available
270 """
271 ref=-Lsup(f)*2
272 f_max=ref
273 t_max=None
274 p_max=None
275 x=f.getFunctionSpace().getX()
276 for t in self.getTags():
277 p,m=self.getParametrization(x,tag=t, tol=tol)
278 loc=((m*f)+(1.-m)*ref).maxGlobalDataPoint()
279 f_t=f.getTupleForGlobalDataPoint(*loc)[0]
280 if f_t>f_max:
281 f_max=f_t
282 t_max=t
283 p_max=p.getTupleForGlobalDataPoint(*loc)[0]
284
285 return f_max, t_max, p_max
286
287 def getParametrization(self,x,tag=None, tol=sqrt(EPSILON), outsider=None):
288 """
289 returns the parametrization of the fault ``tag`` in the fault system. In fact the values of the parametrization for at given coordinates ``x`` is returned. In addition to the value of the parametrization a mask is returned indicating if the given location is on the fault with given tolerance ``tol``.
290
291 Typical usage of the this method is
292
293 dom=Domain(..)
294 x=dom.getX()
295 fs=FaultSystem()
296 fs.addFault(tag=3,...)
297 p, m=fs.getParametrization(x, outsider=0,tag=3)
298 saveDataCSV('x.csv',p=p, x=x, mask=m)
299
300 to create a file with the coordinates of the points in ``x`` which are on the fault (as ``mask=m``) together with their location ``p`` in the fault coordinate system.
301
302 :param x: location(s)
303 :type x: `escript.Data` object or `numpy.ndarray`
304 :param tag: the tage of the fault
305 :param tol: relative tolerance to check if location is on fault.
306 :type tol: ``float``
307 :param outsider: value used for parametrization values outside the fault. If not present an appropriate value is choosen.
308 :type outsider: ``float``
309 :return: the coordinates ``x`` in the coordinate system of the fault and a mask indicating coordinates in the fault by 1 (0 elsewhere)
310 :rtype: `escript.Data` object or `numpy.ndarray`
311 """
312 offsets=self.getW0Offsets(tag)
313 w0_range=self.getW0Range(tag)[1]-self.getW0Range(tag)[0]
314 if outsider == None:
315 outsider=min(self.getW0Range(tag)[0],self.getW0Range(tag)[1])-abs(w0_range)/sqrt(EPSILON)
316
317 p=x[0]*0 + outsider
318 updated=x[0]*0
319
320 if self.getDim()==2:
321 #
322 #
323 s=self.getFaultSegments(tag)[0]
324 for i in xrange(1,len(s)):
325 d=s[i]-s[i-1]
326 h=x-s[i-1]
327 h_l=length(h)
328 d_l=length(d)
329 if not d_l>0:
330 raise ValueError,"segement %s in fault %s has zero length."%(i,tag)
331 t=inner(h,d)/d_l**2
332 t=t*whereNonPositive(t-1.)*whereNonNegative(t)
333 r=length(h-t*d)/maximum(h_l,d_l)
334 m=whereNonPositive(r-tol)*(1.-updated)
335 p=(1.-m)*p+m*(offsets[i-1]+(offsets[i]-offsets[i-1])*t)
336 updated=wherePositive(updated+m)
337 else:
338 raise ValueError,"3D is not supported yet."
339
340 return p, updated
341
342
343 def patchMap(v,v1,v2,v3,v4,x_offset,x_length,depth):
344 """
345 maps a patch with corners v1,v2,v3,v4 onto
346 the rectangle [x_offset,0] x [x_offset+x_length,-depth]
347 """
348 #
349 # inspect the mapping v-> w=(w0,w1) with 0<=s<=1 0<=t<=1
350 # then we apply x_offset+s*x_length, -depth*t to map to the final domain
351 #
352 # we need to find s and t such that
353 #
354 # v=v1*(1-t)*(1-s)+v2*t*(1-s)+v3*t*s+v4*(1-t)*s
355 # =v1*(1-t-s+t*s)+v2*(t-t*s)+v3*t*s+v4*(s-t*s)
356 # =v1+t*(v2-v1)+s*(v4-v1)+s*t*(v1-v2+v3-v4)
357 #
358 # with d4=(v2-v1)-dot(v2-v1,v4-v1)/dot(v4-v1,v4-v1)*(v4-v1) => dot(d4,v2-v1)=0
359 # with d2=(v4-v1)-dot(v2-v1,v4-v1)/dot(v2-v1,v2-v1)*(v2-v1) => dot(d2,v4-v1)=0
360 #
361 # dot(v-v1,d4) = s*dot(v4-v1,d4)+s*t*dot(v3-v4,d4) (a)
362 # dot(v-v1,d2) = t*dot(v2-v1,d2)+s*t*dot(v3-v2,d2) (b)
363 #
364 # in addition we set n=cross(d2,d4) => dot(v2-v1,n)=dot(v4-v1,n)=0
365 #
366 # dot(v-v1,n)=s*t*dot(v3-v1,n)
367 #
368 # => s*t=dot(v-v1,n)/dot(v3-v1,n) if dot(v3-v1,n)!=0
369 # otherwise on can set s*t=0 to solve (a) and (b)
370 #
371 h=v-v1
372 h2=(v2-v1)
373 h3=(v3-v1)
374 h4=(v4-v1)
375 d4=h4-numpy.dot(h2,h4)/numpy.dot(h2,h2)*h2
376 d2=h2-numpy.dot(h2,h4)/numpy.dot(h4,h4)*h4
377 n=numpy.cross(d4,d2)
378 h31n=numpy.dot(n,v3-v1)
379 if abs(h31n)>0:
380 st=n/h31n
381 else:
382 st=0*n
383 s=numpy.dot(h,d4-st*numpy.dot(v3-v4,d4))/numpy.dot(v4-v1,d4)
384 t=numpy.dot(h,d2-st*numpy.dot(v3-v2,d2))/numpy.dot(v2-v1,d2)
385 return x_offset+s*x_length, -depth*t

  ViewVC Help
Powered by ViewVC 1.1.26