/[escript]/trunk/esys2/escript/py_src/util.py
ViewVC logotype

Contents of /trunk/esys2/escript/py_src/util.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 122 - (show annotations)
Thu Jun 9 05:38:05 2005 UTC (14 years, 2 months ago) by jgs
File MIME type: text/x-python
File size: 12104 byte(s)
Merge of development branch back to main trunk on 2005-06-09

1 # $Id$
2
3 ## @file util.py
4
5 """
6 @brief Utility functions for escript
7 """
8
9 import numarray
10 import escript
11 #
12 # escript constants (have to be consistent witj utilC.h
13 #
14 UNKNOWN=-1
15 EPSILON=1.e-15
16 Pi=numarray.pi
17 # some solver options:
18 NO_REORDERING=0
19 MINIMUM_FILL_IN=1
20 NESTED_DISSECTION=2
21 # solver methods
22 DEFAULT_METHOD=0
23 DIRECT=1
24 CHOLEVSKY=2
25 PCG=3
26 CR=4
27 CGS=5
28 BICGSTAB=6
29 SSOR=7
30 ILU0=8
31 ILUT=9
32 JACOBI=10
33 GMRES=11
34 PRES20=12
35
36 METHOD_KEY="method"
37 SYMMETRY_KEY="symmetric"
38 TOLERANCE_KEY="tolerance"
39
40 # supported file formats:
41 VRML=1
42 PNG=2
43 JPEG=3
44 JPG=3
45 PS=4
46 OOGL=5
47 BMP=6
48 TIFF=7
49 OPENINVENTOR=8
50 RENDERMAN=9
51 PNM=10
52
53 #
54 # wrapper for various functions: if the argument has attribute the function name
55 # as an argument it calls the correspong methods. Otherwise the coresponsing numarray
56 # function is called.
57 #
58 # functions involving the underlying Domain:
59 #
60 def grad(arg,where=None):
61 """
62 @brief returns the spatial gradient of the Data object arg
63
64 @param arg: Data object representing the function which gradient to be calculated.
65 @param where: FunctionSpace in which the gradient will be. If None Function(dom) where dom is the
66 domain of the Data object arg.
67 """
68 if isinstance(arg,escript.Data):
69 if where==None:
70 return arg.grad()
71 else:
72 return arg.grad(where)
73 else:
74 return arg*0.
75
76 def integrate(arg,what=None):
77 """
78 @brief return the integral if the function represented by Data object arg over its domain.
79
80 @param arg
81 """
82 if not what==None:
83 arg2=escript.Data(arg,what)
84 else:
85 arg2=arg
86 if arg2.getRank()==0:
87 return arg2.integrate()[0]
88 else:
89 return arg2.integrate()
90
91 def interpolate(arg,where):
92 """
93 @brief interpolates the function represented by Data object arg into the FunctionSpace where.
94
95 @param arg
96 @param where
97 """
98 if isinstance(arg,escript.Data):
99 return arg.interpolate(where)
100 else:
101 return arg
102
103 # functions returning Data objects:
104
105 def transpose(arg,axis=None):
106 """
107 @brief returns the transpose of the Data object arg.
108
109 @param arg
110 """
111 if isinstance(arg,escript.Data):
112 # hack for transpose
113 r=arg.getRank()
114 if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"
115 s=arg.getShape()
116 out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())
117 for i in range(s[0]):
118 for j in range(s[1]):
119 out[j,i]=arg[i,j]
120 return out
121 # end hack for transpose
122 if axis==None: axis=arg.getRank()/2
123 return arg.transpose(axis)
124 else:
125 if axis==None: axis=arg.rank/2
126 return numarray.transpose(arg,axis=axis)
127
128 def trace(arg):
129 """
130 @brief
131
132 @param arg
133 """
134 if isinstance(arg,escript.Data):
135 # hack for trace
136 r=arg.getRank()
137 if r!=2: raise ValueError,"trace only avalaible for rank 2 objects"
138 s=arg.getShape()
139 out=escript.Scalar(0,arg.getFunctionSpace())
140 for i in range(min(s)):
141 out+=arg[i,i]
142 return out
143 # end hack for trace
144 return arg.trace()
145 else:
146 return numarray.trace(arg)
147
148 def exp(arg):
149 """
150 @brief
151
152 @param arg
153 """
154 if isinstance(arg,escript.Data):
155 return arg.exp()
156 else:
157 return numarray.exp(arg)
158
159 def sqrt(arg):
160 """
161 @brief
162
163 @param arg
164 """
165 if isinstance(arg,escript.Data):
166 return arg.sqrt()
167 else:
168 return numarray.sqrt(arg)
169
170 def sin(arg):
171 """
172 @brief
173
174 @param arg
175 """
176 if isinstance(arg,escript.Data):
177 return arg.sin()
178 else:
179 return numarray.sin(arg)
180
181 def tan(arg):
182 """
183 @brief
184
185 @param arg
186 """
187 if isinstance(arg,escript.Data):
188 return arg.tan()
189 else:
190 return numarray.tan(arg)
191
192 def cos(arg):
193 """
194 @brief
195
196 @param arg
197 """
198 if isinstance(arg,escript.Data):
199 return arg.cos()
200 else:
201 return numarray.cos(arg)
202
203 def maxval(arg):
204 """
205 @brief
206
207 @param arg
208 """
209 if isinstance(arg,escript.Data):
210 return arg.maxval()
211 elif isinstance(arg,float) or isinstance(arg,int):
212 return arg
213 else:
214 return arg.max()
215
216 def minval(arg):
217 """
218 @brief
219
220 @param arg
221 """
222 if isinstance(arg,escript.Data):
223 return arg.minval()
224 elif isinstance(arg,float) or isinstance(arg,int):
225 return arg
226 else:
227 return arg.min()
228
229 def length(arg):
230 """
231 @brief
232
233 @param arg
234 """
235 if isinstance(arg,escript.Data):
236 if arg.isEmpty(): return escript.Data()
237 if arg.getRank()==0:
238 return abs(arg)
239 elif arg.getRank()==1:
240 sum=escript.Scalar(0,arg.getFunctionSpace())
241 for i in range(arg.getShape()[0]):
242 sum+=arg[i]**2
243 return sqrt(sum)
244 elif arg.getRank()==2:
245 sum=escript.Scalar(0,arg.getFunctionSpace())
246 for i in range(arg.getShape()[0]):
247 for j in range(arg.getShape()[1]):
248 sum+=arg[i,j]**2
249 return sqrt(sum)
250 elif arg.getRank()==3:
251 sum=escript.Scalar(0,arg.getFunctionSpace())
252 for i in range(arg.getShape()[0]):
253 for j in range(arg.getShape()[1]):
254 for k in range(arg.getShape()[2]):
255 sum+=arg[i,j,k]**2
256 return sqrt(sum)
257 elif arg.getRank()==4:
258 sum=escript.Scalar(0,arg.getFunctionSpace())
259 for i in range(arg.getShape()[0]):
260 for j in range(arg.getShape()[1]):
261 for k in range(arg.getShape()[2]):
262 for l in range(arg.getShape()[3]):
263 sum+=arg[i,j,k,l]**2
264 return sqrt(sum)
265 else:
266 raise SystemError,"length is not been implemented yet"
267 # return arg.length()
268 else:
269 return sqrt((arg**2).sum())
270
271 def deviator(arg):
272 """
273 @brief
274
275 @param arg1
276 """
277 if isinstance(arg,escript.Data):
278 shape=arg.getShape()
279 else:
280 shape=arg.shape
281 if len(shape)!=2:
282 raise ValueError,"Deviator requires rank 2 object"
283 if shape[0]!=shape[1]:
284 raise ValueError,"Deviator requires a square matrix"
285 return arg-1./(shape[0]*1.)*trace(arg)*kronecker(shape[0])
286
287 def inner(arg1,arg2):
288 """
289 @brief
290
291 @param arg1, arg2
292 """
293 sum=escript.Scalar(0,arg1.getFunctionSpace())
294 if arg.getRank()==0:
295 return arg1*arg2
296 elif arg.getRank()==1:
297 sum=escript.Scalar(0,arg.getFunctionSpace())
298 for i in range(arg.getShape()[0]):
299 sum+=arg1[i]*arg2[i]
300 elif arg.getRank()==2:
301 sum=escript.Scalar(0,arg.getFunctionSpace())
302 for i in range(arg.getShape()[0]):
303 for j in range(arg.getShape()[1]):
304 sum+=arg1[i,j]*arg2[i,j]
305 elif arg.getRank()==3:
306 sum=escript.Scalar(0,arg.getFunctionSpace())
307 for i in range(arg.getShape()[0]):
308 for j in range(arg.getShape()[1]):
309 for k in range(arg.getShape()[2]):
310 sum+=arg1[i,j,k]*arg2[i,j,k]
311 elif arg.getRank()==4:
312 sum=escript.Scalar(0,arg.getFunctionSpace())
313 for i in range(arg.getShape()[0]):
314 for j in range(arg.getShape()[1]):
315 for k in range(arg.getShape()[2]):
316 for l in range(arg.getShape()[3]):
317 sum+=arg1[i,j,k,l]*arg2[i,j,k,l]
318 else:
319 raise SystemError,"inner is not been implemented yet"
320 return sum
321
322 def sign(arg):
323 """
324 @brief
325
326 @param arg
327 """
328 if isinstance(arg,escript.Data):
329 return arg.sign()
330 else:
331 return numarray.greater(arg,numarray.zeros(arg.shape))-numarray.less(arg,numarray.zeros(arg.shape))
332
333 # reduction operations:
334
335 def sum(arg):
336 """
337 @brief
338
339 @param arg
340 """
341 return arg.sum()
342
343 def sup(arg):
344 """
345 @brief
346
347 @param arg
348 """
349 if isinstance(arg,escript.Data):
350 return arg.sup()
351 elif isinstance(arg,float) or isinstance(arg,int):
352 return arg
353 else:
354 return arg.max()
355
356 def inf(arg):
357 """
358 @brief
359
360 @param arg
361 """
362 if isinstance(arg,escript.Data):
363 return arg.inf()
364 elif isinstance(arg,float) or isinstance(arg,int):
365 return arg
366 else:
367 return arg.min()
368
369 def L2(arg):
370 """
371 @brief returns the L2-norm of the
372
373 @param arg
374 """
375 if isinstance(arg,escript.Data):
376 return arg.L2()
377 elif isinstance(arg,float) or isinstance(arg,int):
378 return abs(arg)
379 else:
380 return numarry.sqrt(dot(arg,arg))
381
382 def Lsup(arg):
383 """
384 @brief
385
386 @param arg
387 """
388 if isinstance(arg,escript.Data):
389 return arg.Lsup()
390 elif isinstance(arg,float) or isinstance(arg,int):
391 return abs(arg)
392 else:
393 return max(numarray.abs(arg))
394
395 def Linf(arg):
396 """
397 @brief
398
399 @param arg
400 """
401 if isinstance(arg,escript.Data):
402 return arg.Linf()
403 elif isinstance(arg,float) or isinstance(arg,int):
404 return abs(arg)
405 else:
406 return min(numarray.abs(arg))
407
408 def dot(arg1,arg2):
409 """
410 @brief
411
412 @param arg
413 """
414 if isinstance(arg1,escript.Data):
415 return arg1.dot(arg2)
416 elif isinstance(arg1,escript.Data):
417 return arg2.dot(arg1)
418 else:
419 return numarray.dot(arg1,arg2)
420
421 def kronecker(d):
422 if hasattr(d,"getDim"):
423 return numarray.identity(d.getDim())
424 else:
425 return numarray.identity(d)
426
427 def unit(i,d):
428 """
429 @brief return a unit vector of dimension d with nonzero index i
430 @param d dimension
431 @param i index
432 """
433 e = numarray.zeros((d,))
434 e[i] = 1.0
435 return e
436
437 #
438 # $Log$
439 # Revision 1.10 2005/06/09 05:37:59 jgs
440 # Merge of development branch back to main trunk on 2005-06-09
441 #
442 # Revision 1.2.2.14 2005/05/20 04:05:23 gross
443 # some work on a darcy flow started
444 #
445 # Revision 1.2.2.13 2005/03/16 05:17:58 matt
446 # Implemented unit(idx, dim) to create cartesian unit basis vectors to
447 # complement kronecker(dim) function.
448 #
449 # Revision 1.2.2.12 2005/03/10 08:14:37 matt
450 # Added non-member Linf utility function to complement Data::Linf().
451 #
452 # Revision 1.2.2.11 2005/02/17 05:53:25 gross
453 # some bug in saveDX fixed: in fact the bug was in
454 # DataC/getDataPointShape
455 #
456 # Revision 1.2.2.10 2005/01/11 04:59:36 gross
457 # automatic interpolation in integrate switched off
458 #
459 # Revision 1.2.2.9 2005/01/11 03:38:13 gross
460 # Bug in Data.integrate() fixed for the case of rank 0. The problem is not totallly resolved as the method should return a scalar rather than a numarray object in the case of rank 0. This problem is fixed by the util.integrate wrapper.
461 #
462 # Revision 1.2.2.8 2005/01/05 04:21:41 gross
463 # FunctionSpace checking/matchig in slicing added
464 #
465 # Revision 1.2.2.7 2004/12/29 05:29:59 gross
466 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
467 #
468 # Revision 1.2.2.6 2004/12/24 06:05:41 gross
469 # some changes in linearPDEs to add AdevectivePDE
470 #
471 # Revision 1.2.2.5 2004/12/17 00:06:53 gross
472 # mk sets ESYS_ROOT is undefined
473 #
474 # Revision 1.2.2.4 2004/12/07 03:19:51 gross
475 # options for GMRES and PRES20 added
476 #
477 # Revision 1.2.2.3 2004/12/06 04:55:18 gross
478 # function wraper extended
479 #
480 # Revision 1.2.2.2 2004/11/22 05:44:07 gross
481 # a few more unitary functions have been added but not implemented in Data yet
482 #
483 # Revision 1.2.2.1 2004/11/12 06:58:15 gross
484 # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
485 #
486 # Revision 1.2 2004/10/27 00:23:36 jgs
487 # fixed minor syntax error
488 #
489 # Revision 1.1.1.1 2004/10/26 06:53:56 jgs
490 # initial import of project esys2
491 #
492 # Revision 1.1.2.3 2004/10/26 06:43:48 jgs
493 # committing Lutz's and Paul's changes to brach jgs
494 #
495 # Revision 1.1.4.1 2004/10/20 05:32:51 cochrane
496 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
497 #
498 # Revision 1.1 2004/08/05 03:58:27 gross
499 # Bug in Assemble_NodeCoordinates fixed
500 #
501 #

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26