35 |
__version__="$Revision:$" |
__version__="$Revision:$" |
36 |
__date__="$Date:$" |
__date__="$Date:$" |
37 |
|
|
38 |
from esys.escript import Lsup,whereZero,kronecker |
from esys.escript import * |
39 |
from esys.escript.benchmark import BenchmarkProblem, Options, BenchmarkFilter |
from esys.escript.benchmark import BenchmarkProblem, Options, BenchmarkFilter |
40 |
import esys.finley |
import esys.finley |
41 |
from esys.escript.linearPDEs import LinearPDE |
from esys.escript.linearPDEs import LinearPDE |
42 |
import os |
import os |
43 |
import math |
import math |
44 |
|
import numarray |
45 |
|
|
46 |
class FinleyFilter(BenchmarkFilter): |
class FinleyFilter(BenchmarkFilter): |
47 |
""" |
""" |
313 |
""" |
""" |
314 |
def __init__(self,n,order,dim,mu0,normal,alpha): |
def __init__(self,n,order,dim,mu0,normal,alpha): |
315 |
self.mu0=mu0 |
self.mu0=mu0 |
316 |
self.normal=normal |
self.normal=numarray.array(normal) |
317 |
self.alpha=alpha |
self.alpha=alpha |
318 |
super(AnisotropicSystem,self).__init__(n,order,dim) |
super(AnisotropicSystem,self).__init__(n,order,dim) |
319 |
|
|
333 |
msk=whereZero(x[0])+whereZero(x[0]-1.) |
msk=whereZero(x[0])+whereZero(x[0]-1.) |
334 |
for i in range(1,d): |
for i in range(1,d): |
335 |
msk+=whereZero(x[i])+whereZero(x[i]-1.) |
msk+=whereZero(x[i])+whereZero(x[i]-1.) |
336 |
msk=msk*numarry.ones((d,),numarray.float) |
msk=msk*numarray.ones((d,),numarray.Float) |
337 |
|
|
338 |
u=x[:] |
u=x[:] |
339 |
for i in range(d): |
for i in range(d): |
342 |
if not i==k: s=s*x[k] |
if not i==k: s=s*x[k] |
343 |
u[i]+=1./d*s |
u[i]+=1./d*s |
344 |
|
|
345 |
s=whereNegative(inner(x-numarry.ones((d,),numarray.float)/2,self.normal)) |
s=whereNegative(inner(x-numarray.ones((d,),numarray.Float)/2,self.normal)) |
346 |
mu=s+mu0*(1.-s) |
mu=s+self.mu0*(1.-s) |
347 |
lam=max(1.,self.mu0)*self.alpha |
lam=max(1.,self.mu0)*self.alpha |
348 |
|
|
349 |
F=Tensor(0.,Function(domain)) |
F=Tensor(0.,Function(domain)) |
358 |
s*=x[k] |
s*=x[k] |
359 |
F[i,j]+=mu/d*s |
F[i,j]+=mu/d*s |
360 |
C=Tensor4(0.,Function(domain)) |
C=Tensor4(0.,Function(domain)) |
361 |
for i in range(w.getDim()): |
for i in range(domain.getDim()): |
362 |
for j in range(w.getDim()): |
for j in range(domain.getDim()): |
363 |
C[i,i,j,j]+=lam |
C[i,i,j,j]+=lam |
364 |
C[j,i,j,i]+=mu |
C[j,i,j,i]+=mu |
365 |
C[j,i,i,j]+=mu |
C[j,i,i,j]+=mu |