source: inundation/ga/storm_surge/pyvolution/util.py @ 195

Last change on this file since 195 was 195, checked in by ole, 20 years ago
File size: 9.2 KB
Line 
1
2def distance(v0, v1):
3    from math import sqrt
4
5    return sqrt((v1[0]-v0[0])**2 + (v1[1]-v0[1])**2)
6
7def length(v):
8    from math import sqrt
9
10    return sqrt(v[0]**2 + v[1]**2)
11
12def inner_product(v0, v1):
13    from math import sqrt
14
15    #FIXME: Obsolete
16    try:
17        i = sqrt(v0[0]*v1[0] + v0[1]*v1[1])
18    except:
19        print v0, v1
20        import sys; sys.exit()
21
22    return i   
23
24def angle(v):
25    """Compute angle between e1 (the unit vector in the x-direction)
26    and the specified vector
27    """
28
29    from math import acos, pi
30
31    l = length(v)
32    v1 = v[0]/l
33    v2 = v[1]/l
34       
35    theta = acos(v1)
36
37    if v2 < 0:
38        #Quadrant 3 or 4
39        theta = 2*pi-theta
40
41    return theta
42
43
44def anglediff(v0, v1):
45    """Compute difference between angle of vector x0, y0 and x1, y1.
46    This is used for determining the ordering of vertices,
47    e.g. for checking if they are counter clockwise.
48   
49    Always return a positive value
50    """
51       
52    from math import pi
53   
54    a0 = angle(v0)
55    a1 = angle(v1)
56
57    #Ensure that difference will be positive
58    if a0 < a1:
59        a0 += 2*pi
60           
61    return a0-a1
62
63
64 
65def compile(FNs=None, CC=None, LD = None, SFLAG = None, verbose = 1):
66  """compile(FNs=None, CC=None, LD = None, SFLAG = None):
67 
68     Compile FN(s) using compiler CC (e.g. mpicc),
69     Loader LD and shared flag SFLAG.
70     If CC is absent use default compiler dependent on platform
71     if LD is absent CC is used.
72     if SFLAG is absent platform default is used
73     FNs can be either one filename or a list of filenames
74     In the latter case, the first will be used to name so file.
75  """
76  import os, string, sys, types
77 
78  # Input check
79  #
80  assert not FNs is None, "No filename provided"
81
82  if not type(FNs) == types.ListType:
83    FNs = [FNs]
84
85
86  libext = 'so' #Default extension (Unix)
87  libs = ''
88  version = sys.version[:3]
89 
90  # Determine platform and compiler
91  #
92  if sys.platform == 'sunos5':  #Solaris
93    if CC:
94      compiler = CC
95    else: 
96      compiler = 'gcc'
97    if LD:
98      loader = LD
99    else: 
100      loader = compiler
101    if SFLAG:
102      sharedflag = SFLAG
103    else: 
104      sharedflag = 'G'
105     
106  elif sys.platform == 'osf1V5':  #Compaq AlphaServer
107    if CC:
108      compiler = CC
109    else: 
110      compiler = 'cc'
111    if LD:
112      loader = LD
113    else: 
114      loader = compiler
115    if SFLAG:
116      sharedflag = SFLAG
117    else: 
118      sharedflag = 'shared'   
119     
120  elif sys.platform == 'linux2':  #Linux
121    if CC:
122      compiler = CC
123    else: 
124      compiler = 'gcc'
125    if LD:
126      loader = LD
127    else: 
128      loader = compiler
129    if SFLAG:
130      sharedflag = SFLAG
131    else: 
132      sharedflag = 'shared'   
133     
134  elif sys.platform == 'darwin':  #Mac OS X:
135    if CC:
136      compiler = CC
137    else: 
138      compiler = 'cc'
139    if LD:
140      loader = LD
141    else: 
142      loader = compiler
143    if SFLAG:
144      sharedflag = SFLAG
145    else: 
146      sharedflag = 'bundle -flat_namespace -undefined suppress'
147
148  elif sys.platform == 'cygwin':  #Cygwin (compilation same as linux)
149    if CC:
150      compiler = CC
151    else: 
152      compiler = 'gcc'
153    if LD:
154      loader = LD
155    else: 
156      loader = compiler
157    if SFLAG:
158      sharedflag = SFLAG
159    else: 
160      sharedflag = 'shared'
161    libext = 'dll'
162    libs = '/lib/python%s/config/libpython%s.dll.a' %(version,version)
163     
164  elif sys.platform == 'win32':  #Windows
165    if CC:
166      compiler = CC
167    else: 
168      compiler = 'gcc'
169    if LD:
170      loader = LD
171    else: 
172      loader = compiler
173    if SFLAG:
174      sharedflag = SFLAG
175    else: 
176      sharedflag = 'shared'
177    libext = 'dll'
178
179    v = version.replace('.','')
180    dllfilename = 'python%s.dll' %(v)
181    libs = os.path.join(sys.exec_prefix,dllfilename)
182     
183     
184  else:
185    if verbose: print "Unrecognised platform %s - revert to default"\
186                %sys.platform
187   
188    if CC:
189      compiler = CC
190    else: 
191      compiler = 'cc'
192    if LD:
193      loader = LD
194    else: 
195      loader = 'ld'
196    if SFLAG:
197      sharedflag = SFLAG
198    else: 
199      sharedflag = 'G'
200
201           
202  # Find location of include files
203  #
204  if sys.platform == 'win32':  #Windows
205    include = os.path.join(sys.exec_prefix, 'include')   
206  else: 
207    include = os.path.join(os.path.join(sys.exec_prefix, 'include'),
208                           'python'+version)
209
210  # Check existence of Python.h
211  #
212  headerfile = include + os.sep +'Python.h'
213  try:
214    open(headerfile, 'r')
215  except:
216    raise """Did not find Python header file %s.
217    Make sure files for Python C-extensions are installed.
218    In debian linux, for example, you need to install a
219    package called something like python2.3-dev""" %headerfile
220 
221 
222  # Check filename(s)
223  #
224  object_files = ''
225  for FN in FNs:       
226    root, ext = os.path.splitext(FN)
227    if ext == '':
228      FN = FN + '.c'
229    elif ext.lower() != '.c':
230      raise Exception, "Unrecognised extension: " + FN
231   
232    try:
233      open(FN,'r')
234    except:   
235      raise Exception, "Could not open: " + FN
236
237    if not object_files: root1 = root  #Remember first filename       
238    object_files += root + '.o ' 
239 
240 
241    # Compile
242    #
243   
244    s = "%s -c %s -I%s -o %s.o -Wall -O" %(compiler, FN, include, root)
245    if verbose:
246      print s
247    else:
248      s = s + ' 2> /dev/null' #Suppress errors
249 
250    try:
251      err = os.system(s)
252      if err != 0:
253          raise 'Attempting to compile %s failed - please try manually' %FN
254    except:
255      raise 'Could not compile %s - please try manually' %FN 
256
257 
258  # Make shared library (*.so or *.dll)
259 
260  s = "%s -%s %s -o %s.%s %s -lm" %(loader, sharedflag, object_files, root1, libext, libs)
261  if verbose:
262    print s
263  else:
264    s = s + ' 2> /dev/null' #Suppress warnings
265 
266  try: 
267    err=os.system(s)
268    if err != 0:       
269        raise 'Atempting to link %s failed - please try manually' %root1     
270  except:
271    raise 'Could not link %s - please try manually' %root1
272   
273
274
275
276def can_use_C_extension(filename):
277    """Determine whether specified C-extension
278    can and should be used.
279    """
280
281    from config import use_extensions
282
283    from os.path import splitext
284
285    root, ext = splitext(filename)
286   
287    C=False
288    if use_extensions:
289        try:
290            s = 'import %s' %root
291            #print s
292            exec(s)
293        except:
294            try:
295                open(filename)
296            except:
297                msg = 'C extension %s cannot be opened' %filename
298                print msg               
299            else:   
300                print '------- Trying to compile c-extension %s' %filename
301           
302                compile(filename)                   
303                try:
304                    compile(filename)
305                except:
306                    print 'WARNING: Could not compile C-extension %s'\
307                          %filename
308                else:
309                    try:
310                        exec('import %s' %root)
311                    except:
312                        msg = 'C extension %s seems to compile OK, '
313                        msg += 'but it can still not be imported.'
314                        raise msg
315                    else:
316                        C=True
317        else:
318            C=True
319           
320    if not C:
321        pass
322        print 'NOTICE: C-extension %s not used' %filename
323
324    return C
325
326
327
328
329
330####################################################################
331#Python versions of function that are also implemented in util_ext.c
332#
333
334def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
335    """
336    """
337   
338    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)           
339    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
340    a /= det
341
342    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
343    b /= det           
344
345    return a, b
346
347
348
349def rotate_python(q, normal, direction = 1):
350    """Rotate the momentum component q (q[1], q[2])
351    from x,y coordinates to coordinates based on normal vector.
352
353    If direction is negative the rotation is inverted.
354   
355    Input vector is preserved
356
357    This function is specific to the shallow water wave equation
358    """
359
360    #FIXME: Needs to be tested
361
362    from Numeric import zeros, Float
363   
364    assert len(q) == 3,\
365           'Vector of conserved quantities must have length 3'\
366           'for 2D shallow water equation'
367    try:
368        l = len(normal)
369    except:
370        raise 'Normal vector must be an Numeric array'
371
372    assert l == 2, 'Normal vector must have 2 components'
373
374 
375    n1 = normal[0]
376    n2 = normal[1]   
377   
378    r = zeros(len(q), Float) #Rotated quantities
379    r[0] = q[0]              #First quantity, height, is not rotated
380
381    if direction == -1:
382        n2 = -n2
383
384
385    r[1] =  n1*q[1] + n2*q[2]
386    r[2] = -n2*q[1] + n1*q[2]
387   
388    return r
389
390
391
392
393##############################################
394#Initialise module
395
396import util
397if util.can_use_C_extension('util_ext.c'):
398    from util_ext import gradient, rotate
399else:
400    gradient = gradient_python
401    rotate = rotate_python
402
403
404
405
406if __name__ == "__main__":
407    pass
408
Note: See TracBrowser for help on using the repository browser.