« back to all changes in this revision
Viewing changes to FIAT/expansions.py
- browse files at revision 2
- compare with another revision
- download diff
- download tarball
- view history from revision 2
- Committer: kirby
- Date: 2005-05-31 15:03:57 UTC
- Revision ID: devnull@localhost-20050531150357-n7ka5987yo6ylpaz
[project @ 2005-05-31 15:03:57 by kirby]
Initial revision
Initial revision
- files added:
- FIAT
- doc
added
removed
FIAT/expansions.py
1
# Written by Robert C. Kirby
2
# Copyright 2005 by The University of Chicago
3
# Distributed under the LGPL license
4
# This work is partially supported by the US Department of Energy
5
# under award number DE-FG02-04ER25650
6
7
# Last modified 4 Feb 2005 by RCK
8
9
import jacobi, shapes, math, Numeric
10
11
"""Principal expansion functions as defined by Karniadakis & Sherwin.
12
The main point of entry to this module is the function
13
make_expansions( shape , n ) which takes the simplex shape and the
14
polynomal degree n and returns the list of expansion functions up to
15
that degree."""
16
17
18
def psitilde_a( i , z ):
19
return jacobi.eval_jacobi( 0 , 0 , i , z )
20
21
def psitilde_b( i , j , z ):
22
if i == 0:
23
return jacobi.eval_jacobi( 1 , 0 , j , z )
24
else:
25
return (0.5 * ( 1 - z )) ** i * jacobi.eval_jacobi( 2*i+1 , 0 , j , z )
26
27
def psitilde_c( i , j , k , z ):
28
if i + j == 0:
29
return jacobi.eval_jacobi( 2 , 0 , k , z )
30
else:
31
return ( 0.5 * (1-z))**(i+j) \
32
* jacobi.eval_jacobi( 2*(i+j+1), 0 , k , z )
33
34
# coordinate changes from Karniadakis & Sherwin
35
36
def make_scalings( n , etas ):
37
scalings = Numeric.zeros( (n+1,len(etas)) ,"d")
38
scalings[0,:] = 1.0
39
if n > 0:
40
scalings[1,:] = 0.5 * (1.0 - etas)
41
for k in range(2,n+1):
42
scalings[k,:] = scalings[k-1,:] * scalings[1,:]
43
return scalings
44
45
46
def eta_triangle( xi ):
47
( xi1 , xi2 ) = xi
48
if xi2 == 1.:
49
eta1 = -1.0
50
else:
51
eta1 = 2.0 * ( 1. + xi1 ) / ( 1. - xi2 ) - 1
52
eta2 = xi2
53
return eta1 , eta2
54
55
def xi_triangle( eta ):
56
eta1, eta2 = eta
57
xi1 = 0.5 * (1. + eta1) * (1. - eta2) - 1.
58
xi2 = eta2
59
return xi1,xi2
60
61
def eta_tetrahedron( xi ):
62
xi1,xi2,xi3 = xi
63
if xi2 + xi3 == 0.:
64
eta1 = 1.
65
else:
66
eta1 = -2. * ( 1. + xi1 ) / (xi2 + xi3) - 1.
67
if xi3 == 1.:
68
eta2 = -1.
69
else:
70
eta2 = 2. * (1. + xi2) / (1. - xi3 ) - 1.
71
eta3 = xi3
72
return eta1,eta2,eta3
73
74
75
def xi_tetrahedron( eta ):
76
eta1,eta2,eta3 = eta
77
78
xi1 = 0.25 * ( 1. + eta1 ) * ( 1. - eta2 ) * ( 1. - eta3 ) - 1.
79
xi2 = 0.5 * ( 1. + eta2 ) * ( 1. - eta3 ) - 1.
80
xi3 = eta3
81
82
return xi1,xi2,xi3
83
84
coord_changes = { shapes.TRIANGLE: eta_triangle , \
85
shapes.TETRAHEDRON: eta_tetrahedron }
86
87
inverse_coord_changes = { shapes.TRIANGLE: xi_triangle , \
88
shapes.TETRAHEDRON: xi_tetrahedron }
89
90
def make_coordinate_change( shape ):
91
"""Maps from reference domain to rectangular reference domain."""
92
global coord_changes
93
if coord_changes.has_key( shape ):
94
return coord_changes[shape]
95
else:
96
raise RuntimeError, "Can't collapse coordinates"
97
98
def make_inverse_coordinate_change( shape ):
99
"""Maps from rectangular reference domain to reference domain."""
100
global inverse_coord_changes
101
if inverse_coord_changes.has_key( shape ):
102
return inverse_coord_changes[shape]
103
else:
104
raise RuntimeError, "Can't collapse coordinates"
105
106
def phi_line( index , xi ):
107
(p,) = index
108
(eta,) = xi
109
return psitilde_a( p , eta )
110
111
def phi_triangle( index , xi ):
112
p,q = index
113
eta1,eta2 = eta_triangle( xi )
114
115
alpha = psitilde_a( p , eta1 )
116
beta = psitilde_b( p , q , eta2 )
117
118
return alpha * beta
119
120
def phi_tetrahedron( index , xi ):
121
p,q,r = index
122
eta1,eta2,eta3 = eta_tetrahedron( xi )
123
124
alpha = psitilde_a( p , eta1 )
125
beta = psitilde_b( p , q , eta2 )
126
gamma = psitilde_c( p , q , r , eta3 )
127
128
return alpha * beta * gamma
129
130
def tabulate_phis_line( n , xs ):
131
"""Tabulates all the basis functions over the line up to
132
degree in at points xs."""
133
# extract raw points from singletons
134
ys = Numeric.array( [ x for (x,) in xs ] )
135
phis = jacobi.eval_jacobi_batch(0,0,n,ys)
136
for i in range(0,len(phis)):
137
phis[i,:] *= math.sqrt(1.0*i+0.5)
138
return phis
139
140
def tabulate_phis_derivs_line( n , xs ):
141
ys = Numeric.array( [ x for (x,) in xs ] )
142
phi_derivs = jacobi.eval_jacobi_deriv_batch(0,0,n,ys)
143
for i in range(0,len(xs)):
144
phi_derivs[i,:] *= math.sqrt(1.0*i+0.5)
145
return (phi_derivs,)
146
147
def tabulate_phis_triangle( n , xs ):
148
"""Tabulates all the basis functions over the triangle
149
up to degree n"""
150
151
# unpack coordinates
152
etas = map( eta_triangle , xs )
153
eta1s = Numeric.array( [ eta1 for (eta1,eta2) in etas ] )
154
eta2s = Numeric.array( [ eta2 for (eta1,eta2) in etas ] )
155
156
# get Legendre functions for eta1 direction
157
psitilde_as = jacobi.eval_jacobi_batch(0,0,n,eta1s)
158
159
# scalings ( (1-eta2) / 2 ) ** i
160
scalings = make_scalings( n , eta2s )
161
162
# for i == 0, I can have j == 0...degree
163
# for i == 1, I can have j == 0...degree-1
164
165
# phi_{i,j} requires p_j^{2i+1,0}
166
psitilde_bs = [ jacobi.eval_jacobi_batch(2*i+1,0,n-i,eta2s) \
167
for i in range(0,n+1) ]
168
169
results = Numeric.zeros( (shapes.poly_dims[shapes.TRIANGLE](n),len(xs)) , \
170
"d" )
171
cur = 0
172
for k in range(0,n+1):
173
for i in range(0,k+1):
174
ii = k-i
175
jj = i
176
results[cur,:] = psitilde_as[ii,:] * scalings[ii,:] \
177
* psitilde_bs[ii][jj,:]
178
results[cur,:] *= math.sqrt( (ii+0.5)*(ii+jj+1.0) )
179
cur += 1
180
181
return results
182
183
def tabulate_phis_derivs_triangle(n,xs):
184
"""I'm not properly handling the singularity at the top,
185
so we can't differentiation for eta2 == 0."""
186
# unpack coordinates
187
etas = map( eta_triangle , xs )
188
eta1s = Numeric.array( [ eta1 for (eta1,eta2) in etas ] )
189
eta2s = Numeric.array( [ eta2 for (eta1,eta2) in etas ] )
190
191
psitilde_as = jacobi.eval_jacobi_batch(0,0,n,eta1s)
192
psitilde_derivs_as = jacobi.eval_jacobi_deriv_batch(0,0,n,eta1s)
193
psitilde_bs = [ jacobi.eval_jacobi_batch(2*i+1,0,n-i,eta2s) \
194
for i in range(0,n+1) ]
195
psitilde_derivs_bs = [ jacobi.eval_jacobi_deriv_batch(2*i+1,0,n-i,eta2s) \
196
for i in range(0,n+1) ]
197
198
scalings = Numeric.zeros( (n+1,len(xs)) ,"d")
199
scalings[0,:] = 1.0
200
if n > 0:
201
scalings[1,:] = 0.5 * (1.0 - eta2s)
202
for k in range(2,n+1):
203
scalings[k,:] = scalings[k-1,:] * scalings[1,:]
204
205
xderivs = Numeric.zeros( (shapes.poly_dims[shapes.TRIANGLE](n),len(xs)) , \
206
"d" )
207
yderivs = Numeric.zeros( xderivs.shape , "d" )
208
tmp = Numeric.zeros( (len(xs),) , "d" )
209
210
cur = 0
211
for k in range(0,n+1):
212
for i in range(0,k+1):
213
ii = k-i
214
jj = i
215
216
xderivs[cur,:] = psitilde_derivs_as[ii,:] \
217
* psitilde_bs[ii][jj,:]
218
if ii > 0:
219
xderivs[cur,:] *= scalings[ii-1,:]
220
221
yderivs[cur,:] = psitilde_derivs_as[ii,:] \
222
* psitilde_bs[ii][jj,:] \
223
* 0.5 * (1.0+eta1s[:])
224
if ii > 0:
225
yderivs[cur,:] *= scalings[ii-1,:]
226
227
228
tmp[:] = psitilde_derivs_bs[ii][jj,:] \
229
* scalings[ii,:]
230
if ii > 0:
231
tmp[:] -= 0.5 * ii * psitilde_bs[ii][jj,:] \
232
* scalings[ii-1]
233
234
yderivs[cur,:] += psitilde_as[ii,:] * tmp
235
236
alpha = math.sqrt( (ii+0.5)*(ii+jj+1.0) )
237
xderivs[cur,:] *= alpha
238
yderivs[cur,:] *= alpha
239
cur += 1
240
241
return (xderivs,yderivs)
242
243
def tabulate_phis_tetrahedron( n , xs ):
244
"""Tabulates all the basis functions over the tetrahedron
245
up to degree n"""
246
# unpack coordinates
247
etas = map( eta_tetrahedron , xs )
248
eta1s = Numeric.array( [ eta1 for (eta1,eta2,eta3) in etas ] )
249
eta2s = Numeric.array( [ eta2 for (eta1,eta2,eta3) in etas ] )
250
eta3s = Numeric.array( [ eta3 for (eta1,eta2,eta3) in etas ] )
251
252
# get Legendre functions for eta1 direction
253
psitilde_as = jacobi.eval_jacobi_batch(0,0,n,eta1s)
254
255
eta2_scalings = Numeric.zeros( (n+1,len(xs)) ,"d")
256
eta2_scalings[0,:] = 1.0
257
if n > 0:
258
eta2_scalings[1,:] = 0.5 * (1.0 - eta2s)
259
for k in range(2,n+1):
260
eta2_scalings[k,:] = eta2_scalings[k-1,:] * eta2_scalings[1,:]
261
262
psitilde_bs = [ jacobi.eval_jacobi_batch(2*i+1,0,n-i,eta2s) \
263
for i in range(0,n+1) ]
264
265
# (0.5*(1-z))**i+j, since for k=0, we can have i+j up to n,
266
# we need same structure as eta2_scalings
267
eta3_scalings = Numeric.zeros( (n+1,len(xs)) ,"d" )
268
eta3_scalings[0,:] = 1.0
269
if n > 0:
270
eta3_scalings[1,:] = 0.5 * (1.0 - eta3s)
271
for k in range(2,n+1):
272
eta3_scalings[k,:] = eta3_scalings[k-1,:] * eta3_scalings[1,:]
273
274
# I need psitilde_c[i][j][k]
275
# for 0<=i+j+k<=n
276
277
psitilde_cs = [ [ jacobi.eval_jacobi_batch(2*(i+j+1),0,\
278
n-i-j,eta3s) \
279
for j in range(0,n+1-i) ] for i in range(0,n+1) ]
280
281
results = Numeric.zeros( (shapes.poly_dims[shapes.TETRAHEDRON](n), \
282
len(xs)) , \
283
"d" )
284
cur = 0
285
for k in range(0,n+1): # loop over degree
286
for i in range(0,k+1):
287
for j in range(0,k-i+1):
288
ii = k-i-j
289
jj = j
290
kk = i
291
results[cur,:] = psitilde_as[ii,:] * \
292
eta2_scalings[ii,:] \
293
* psitilde_bs[ii][jj,:] \
294
* eta3_scalings[ii+jj,:] \
295
* psitilde_cs[ii][jj][kk,:]
296
results[cur,:] *= math.sqrt( (ii+0.5) \
297
* (ii+jj+1.0) \
298
* (ii+jj+kk+1.5) )
299
cur += 1
300
301
return results
302
303
def tabulate_phis_derivs_tetrahedron(n,xs):
304
"""Tabulates all the derivatives of basis functions over the tetrahedron
305
up to degree n at points xs"""
306
# unpack coordinates
307
etas = map( eta_tetrahedron , xs )
308
eta1s = Numeric.array( [ eta1 for (eta1,eta2,eta3) in etas ] )
309
eta2s = Numeric.array( [ eta2 for (eta1,eta2,eta3) in etas ] )
310
eta3s = Numeric.array( [ eta3 for (eta1,eta2,eta3) in etas ] )
311
312
psitilde_as = jacobi.eval_jacobi_batch(0,0,n,eta1s)
313
psitilde_as_derivs = jacobi.eval_jacobi_deriv_batch(0,0,n,eta1s)
314
psitilde_bs = [ jacobi.eval_jacobi_batch(2*i+1,0,n-i,eta2s) \
315
for i in range(0,n+1) ]
316
psitilde_bs_derivs = [ jacobi.eval_jacobi_deriv_batch(2*i+1,0,n-i,eta2s) \
317
for i in range(0,n+1) ]
318
psitilde_cs = [ [ jacobi.eval_jacobi_batch(2*(i+j+1),0,\
319
n-i-j,eta3s) \
320
for j in range(0,n+1-i) ] for i in range(0,n+1) ]
321
psitilde_cs_derivs = [ [ jacobi.eval_jacobi_deriv_batch(2*(i+j+1),0,\
322
n-i-j,eta3s) \
323
for j in range(0,n+1-i) ] for i in range(0,n+1) ]
324
325
eta2_scalings = Numeric.zeros( (n+1,len(xs)) ,"d")
326
eta2_scalings[0,:] = 1.0
327
if n > 0:
328
eta2_scalings[1,:] = 0.5 * (1.0 - eta2s)
329
for k in range(2,n+1):
330
eta2_scalings[k,:] = eta2_scalings[k-1,:] * eta2_scalings[1,:]
331
332
eta3_scalings = Numeric.zeros( (n+1,len(xs)) ,"d" )
333
eta3_scalings[0,:] = 1.0
334
if n > 0:
335
eta3_scalings[1,:] = 0.5 * (1.0 - eta3s)
336
for k in range(2,n+1):
337
eta3_scalings[k,:] = eta3_scalings[k-1,:] \
338
* eta3_scalings[1,:]
339
340
tmp = Numeric.zeros( (len(xs),) , "d" )
341
xderivs = Numeric.zeros( (shapes.poly_dims[shapes.TETRAHEDRON](n), \
342
len(xs)) , \
343
"d" )
344
yderivs = Numeric.zeros( xderivs.shape , "d" )
345
zderivs = Numeric.zeros( xderivs.shape , "d" )
346
347
cur = 0
348
for k in range(0,n+1): # loop over degree
349
for i in range(0,k+1):
350
for j in range(0,k-i+1):
351
ii = k-i-j
352
jj = j
353
kk = i
354
355
xderivs[cur,:] = psitilde_as_derivs[ii,:]
356
xderivs[cur,:] *= psitilde_bs[ii][jj,:]
357
xderivs[cur,:] *= psitilde_cs[ii][jj][kk,:]
358
if ii>0:
359
xderivs[cur,:] *= eta2_scalings[ii-1]
360
if ii+jj>0:
361
xderivs[cur,:] *= eta3_scalings[ii+jj-1]
362
363
#d/deta1 with scalings
364
yderivs[cur,:] = psitilde_as_derivs[ii,:]
365
yderivs[cur,:] *= psitilde_bs[ii][jj,:]
366
yderivs[cur,:] *= psitilde_cs[ii][jj][kk,:]
367
yderivs[cur,:] *= 0.5 * (1.0+eta1s[:])
368
if ii>0:
369
yderivs[cur,:] *= eta2_scalings[ii-1]
370
if ii+jj>0:
371
yderivs[cur,:] *= eta3_scalings[ii+jj-1]
372
373
#tmp will hold d/deta2 term
374
tmp[:] = psitilde_bs_derivs[ii][jj,:]
375
tmp[:] *= eta2_scalings[ii,:]
376
if ii>0:
377
tmp[:] -= 0.5*ii*eta2_scalings[ii-1]*psitilde_bs[ii][jj,:]
378
tmp[:] *= psitilde_as[ii,:]
379
tmp[:] *= psitilde_cs[ii][jj][kk,:]
380
if ii+jj>0:
381
tmp[:] *= eta3_scalings[ii+jj-1]
382
yderivs[cur,:] += tmp[:]
383
384
# zderivative
385
#d/deta1 with scalings
386
zderivs[cur,:] = psitilde_as_derivs[ii,:]
387
zderivs[cur,:] *= psitilde_bs[ii][jj,:]
388
zderivs[cur,:] *= psitilde_cs[ii][jj][kk,:]
389
zderivs[cur,:] *= 0.5 * (1.0+eta1s[:])
390
if ii>0:
391
zderivs[cur,:] *= eta2_scalings[ii-1]
392
if ii+jj>0:
393
zderivs[cur,:] *= eta3_scalings[ii+jj-1]
394
395
#tmp will hold d/deta2 term
396
tmp[:] = psitilde_bs_derivs[ii][jj,:]
397
tmp[:] *= eta2_scalings[ii,:]
398
if ii>0:
399
tmp[:] -= 0.5*ii*eta2_scalings[ii-1]*psitilde_bs[ii][jj,:]
400
tmp[:] *= psitilde_as[ii,:]
401
tmp[:] *= psitilde_cs[ii][jj][kk,:]
402
tmp[:] *= 0.5*(1.0+eta2s[:])
403
if ii+jj>0:
404
tmp[:] *= eta3_scalings[ii+jj-1]
405
zderivs[cur,:] += tmp[:]
406
407
#tmp will hold d/deta3 term
408
tmp[:] = psitilde_cs_derivs[ii][jj][kk,:]
409
tmp[:] *= eta3_scalings[ii+jj,:]
410
if ii+jj>0:
411
tmp[:] -= 0.5*(ii+jj)*psitilde_cs[ii][jj][kk,:] \
412
* eta3_scalings[ii+jj-1,:]
413
tmp[:] *= psitilde_as[ii,:]
414
tmp[:] *= psitilde_bs[ii][jj,:]
415
tmp[:] *= eta2_scalings[ii]
416
zderivs[cur,:] += tmp[:]
417
418
xderivs[cur,:] *= math.sqrt( (ii+0.5) \
419
* (ii+jj+1.0) \
420
* (ii+jj+kk+1.5) )
421
422
yderivs[cur,:] *= math.sqrt( (ii+0.5) \
423
* (ii+jj+1.0) \
424
* (ii+jj+kk+1.5) )
425
426
zderivs[cur,:] *= math.sqrt( (ii+0.5) \
427
* (ii+jj+1.0) \
428
* (ii+jj+kk+1.5) )
429
430
431
cur += 1
432
433
return (xderivs,yderivs,zderivs)
434
435
def tabulate_phis_derivs_tetrahedron_old(n,xs):
436
"""Tabulates all the derivatives of basis functions over the tetrahedron
437
up to degree n at points xs"""
438
# unpack coordinates
439
etas = map( eta_tetrahedron , xs )
440
eta1s = Numeric.array( [ eta1 for (eta1,eta2,eta3) in etas ] )
441
eta2s = Numeric.array( [ eta2 for (eta1,eta2,eta3) in etas ] )
442
eta3s = Numeric.array( [ eta3 for (eta1,eta2,eta3) in etas ] )
443
444
psitilde_as = jacobi.eval_jacobi_batch(0,0,n,eta1s)
445
psitilde_as_derivs = jacobi.eval_jacobi_deriv_batch(0,0,n,eta1s)
446
psitilde_bs = [ jacobi.eval_jacobi_batch(2*i+1,0,n-i,eta2s) \
447
for i in range(0,n+1) ]
448
psitilde_bs_derivs = [ jacobi.eval_jacobi_deriv_batch(2*i+1,0,n-i,eta2s) \
449
for i in range(0,n+1) ]
450
psitilde_cs = [ [ jacobi.eval_jacobi_batch(2*(i+j+1),0,\
451
n-i-j,eta3s) \
452
for j in range(0,n+1-i) ] for i in range(0,n+1) ]
453
psitilde_cs_derivs = [ [ jacobi.eval_jacobi_deriv_batch(2*(i+j+1),0,\
454
n-i-j,eta3s) \
455
for j in range(0,n+1-i) ] for i in range(0,n+1) ]
456
457
eta2_scalings = Numeric.zeros( (n+1,len(xs)) ,"d")
458
eta2_scalings[0,:] = 1.0
459
if n > 0:
460
eta2_scalings[1,:] = 0.5 * (1.0 - eta2s)
461
for k in range(2,n+1):
462
eta2_scalings[k,:] = eta2_scalings[k-1,:] * eta2_scalings[1,:]
463
464
eta3_scalings = Numeric.zeros( (n+1,len(xs)) ,"d" )
465
eta3_scalings[0,:] = 1.0
466
if n > 0:
467
eta3_scalings[1,:] = 0.5 * (1.0 - eta3s)
468
for k in range(2,n+1):
469
eta3_scalings[k,:] = eta3_scalings[k-1,:] \
470
* eta3_scalings[1,:]
471
472
tmp = Numeric.zeros( (len(xs),) , "d" )
473
xderivs = Numeric.zeros( (shapes.poly_dims[shapes.TETRAHEDRON](n), \
474
len(xs)) , \
475
"d" )
476
yderivs = Numeric.zeros( xderivs.shape , "d" )
477
zderivs = Numeric.zeros( xderivs.shape , "d" )
478
479
cur = 0
480
for k in range(0,n+1): # loop over degree
481
for i in range(0,k+1):
482
for j in range(0,k-i+1):
483
ii = k-i-j
484
jj = j
485
kk = i
486
487
# xderivative
488
xderivs[cur,:] = psitilde_as_derivs[ii,:] \
489
* psitilde_bs[ii][jj,:] \
490
* psitilde_cs[ii][jj][kk,:]
491
if ii > 0:
492
xderivs[cur,:] *= eta2_scalings[ii-1,:]
493
if ii+jj>0:
494
xderivs[cur,:] *= eta3_scalings[ii+jj-1,:]
495
496
# yderivative
497
# eta1 derivative term
498
yderivs[cur,:] = 0.5 * (1.0+eta1s) * xderivs[cur,:]
499
500
# eta2 derivative term, start with the "internal"
501
# product rule on the eta2 term itself
502
tmp[:] = psitilde_bs_derivs[ii][jj,:] * eta2_scalings[ii,:]
503
if ii>0:
504
tmp[:] -= 0.5 * ii * psitilde_bs[ii][jj,:] \
505
* eta2_scalings[ii-1,:]
506
tmp[:] *= psitilde_as[ii,:] * psitilde_cs[ii][jj][kk,:]
507
if ii+jj>0:
508
tmp[:] *= eta3_scalings[ii+jj-1,:]
509
yderivs[cur,:] += tmp
510
511
# zderivative
512
# eta1 derivative term
513
zderivs[cur,:] = 0.5 * (1.0+eta1s) * xderivs[cur,:]
514
515
# check this term here...
516
# eta2 derivative term
517
zderivs[cur,:] += tmp * 0.5 * (1.0+eta1s)
518
519
# and this one.
520
# eta3 derivative term
521
tmp[:] = psitilde_cs_derivs[ii][jj][kk,:] * eta3_scalings[ii+jj,:]
522
if ii+jj>0:
523
tmp[:] += psitilde_cs[ii][jj][kk,:] * (ii+jj) * eta3_scalings[ii+jj-1,:]
524
tmp[:] *= psitilde_as[ii,:] * psitilde_bs[ii][jj,:] \
525
* eta2_scalings[ii,:]
526
zderivs[cur,:] += tmp
527
528
xderivs[cur,:] *= math.sqrt( (ii+0.5) \
529
* (ii+jj+1.0) \
530
* (ii+jj+kk+1.5) )
531
532
yderivs[cur,:] *= math.sqrt( (ii+0.5) \
533
* (ii+jj+1.0) \
534
* (ii+jj+kk+1.5) )
535
536
zderivs[cur,:] *= math.sqrt( (ii+0.5) \
537
* (ii+jj+1.0) \
538
* (ii+jj+kk+1.5) )
539
540
541
cur += 1
542
543
return (xderivs,yderivs,zderivs)
544
545
546
547
class ExpansionFunction( object ):
548
def __init__( self , indices , phi , alpha ):
549
self.indices, self.phi,self.alpha = indices , phi , alpha
550
return
551
def __call__( self , x ):
552
if len( x ) != len( self.indices ):
553
raise RuntimeError, "Illegal number of coordinates"
554
return self.alpha * self.phi( self.indices , x )
555
556
class PhiLine( ExpansionFunction ):
557
def __init__( self , i ):
558
ExpansionFunction.__init__( self , \
559
(i , ) , \
560
phi_line , \
561
math.sqrt(1.0*i+0.5) )
562
563
class PhiTriangle( ExpansionFunction ):
564
def __init__( self , i , j ):
565
ExpansionFunction.__init__( self , \
566
( i,j ) , \
567
phi_triangle, \
568
math.sqrt( (i+0.5)*(i+j+1.0) ) )
569
return
570
571
class PhiTetrahedron( ExpansionFunction ):
572
def __init__( self , i , j , k ):
573
ExpansionFunction.__init__( self , (i,j,k) , \
574
phi_tetrahedron , \
575
math.sqrt( (i+0.5)*(i+j+1.0)*(i+j+k+1.5)))
576
return
577
578
def make_phis_line( n ):
579
return [ PhiLine( i ) \
580
for i in range(0,n+1) ]
581
582
def make_phis_triangle( n ):
583
return [ PhiTriangle( k - i , i ) \
584
for k in range( 0 , n + 1 ) \
585
for i in range( 0 , k + 1 ) ]
586
587
def make_phis_tetrahedron( n ):
588
return [ PhiTetrahedron( k - i - j , j , i ) \
589
for k in range( 0 , n + 1 ) \
590
for i in range( 0 , k + 1 ) \
591
for j in range( 0 , k - i + 1 ) ]
592
593
make_phis = { shapes.LINE : make_phis_line , \
594
shapes.TRIANGLE : make_phis_triangle , \
595
shapes.TETRAHEDRON : make_phis_tetrahedron }
596
597
def make_expansion( shape , n ):
598
"""Returns the orthogonal expansion basis on a given shape
599
for polynomials of degree n."""
600
global make_phis
601
try:
602
return make_phis[shape]( n )
603
except:
604
raise shapes.ShapeError, "expansions.make_expansion: Illegal shape"
605
606
tabulators = { shapes.LINE : tabulate_phis_line , \
607
shapes.TRIANGLE : tabulate_phis_triangle , \
608
shapes.TETRAHEDRON : tabulate_phis_tetrahedron }
609
610
deriv_tabulators = { shapes.LINE : tabulate_phis_derivs_line , \
611
shapes.TRIANGLE : tabulate_phis_derivs_triangle , \
612
shapes.TETRAHEDRON : tabulate_phis_derivs_tetrahedron }
613
614
615
Loggerhead is a web-based interface for Breezy
Version: 2.0.1