| 1 |
# $Id$ |
| 2 |
|
| 3 |
""" |
| 4 |
provides a some tools related to PDEs currently includes: |
| 5 |
|
| 6 |
Projector - to project a discontinuous |
| 7 |
|
| 8 |
|
| 9 |
""" |
| 10 |
|
| 11 |
from escript import * |
| 12 |
from linearPDEs import LinearPDE |
| 13 |
import numarray |
| 14 |
|
| 15 |
class Projector: |
| 16 |
"""The Projector is a factory which projects a discontiuous function onto a continuous function on the a given domain""" |
| 17 |
def __init__(self, domain, reduce = True, fast=True): |
| 18 |
""" |
| 19 |
@brief Create a continuous function space projector for a domain. |
| 20 |
|
| 21 |
@param domain Domain of the projection. |
| 22 |
@param reduce Flag to reduce projection order (default is True) |
| 23 |
@param fast Flag to use a fast method based on matrix lumping (default is true) |
| 24 |
""" |
| 25 |
self.__pde = LinearPDE(domain) |
| 26 |
self.__pde.setLumping(fast) |
| 27 |
self.__pde.setSymmetryOn() |
| 28 |
self.__pde.setReducedOrderTo(reduce) |
| 29 |
self.__pde.setValue(D = 1.) |
| 30 |
return |
| 31 |
|
| 32 |
def __del__(self): |
| 33 |
return |
| 34 |
|
| 35 |
def __call__(self, input_data): |
| 36 |
""" |
| 37 |
@brief projects input_data onto a continuous function |
| 38 |
|
| 39 |
@param input_data The input_data to be projected. |
| 40 |
""" |
| 41 |
out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain())) |
| 42 |
if input_data.getRank()==0: |
| 43 |
self.__pde.setValue(Y = input_data) |
| 44 |
out=self.__pde.getSolution() |
| 45 |
elif input_data.getRank()==1: |
| 46 |
for i0 in range(input_data.getShape()[0]): |
| 47 |
self.__pde.setValue(Y = input_data[i0]) |
| 48 |
out[i0]=self.__pde.getSolution() |
| 49 |
elif input_data.getRank()==2: |
| 50 |
for i0 in range(input_data.getShape()[0]): |
| 51 |
for i1 in range(input_data.getShape()[1]): |
| 52 |
self.__pde.setValue(Y = input_data[i0,i1]) |
| 53 |
out[i0,i1]=self.__pde.getSolution() |
| 54 |
elif input_data.getRank()==3: |
| 55 |
for i0 in range(input_data.getShape()[0]): |
| 56 |
for i1 in range(input_data.getShape()[1]): |
| 57 |
for i2 in range(input_data.getShape()[2]): |
| 58 |
self.__pde.setValue(Y = input_data[i0,i1,i2]) |
| 59 |
out[i0,i1,i2]=self.__pde.getSolution() |
| 60 |
else: |
| 61 |
for i0 in range(input_data.getShape()[0]): |
| 62 |
for i1 in range(input_data.getShape()[1]): |
| 63 |
for i2 in range(input_data.getShape()[2]): |
| 64 |
for i3 in range(input_data.getShape()[3]): |
| 65 |
self.__pde.setValue(Y = input_data[i0,i1,i2,i3]) |
| 66 |
out[i0,i1,i2,i3]=self.__pde.getSolution() |
| 67 |
return out |
| 68 |
|
| 69 |
|
| 70 |
class Location: |
| 71 |
"""Location provides a factory to access the values of data objects at a given spatial coordinate x. |
| 72 |
In fact, a Location object finds the sample in the set of samples of a given function space or domain which is closest |
| 73 |
to the given point x. """ |
| 74 |
|
| 75 |
def __init__(self,x,where): |
| 76 |
"""initializes a Location to access values in Data objects on the Doamin or FunctionSpace where for the sample point which |
| 77 |
closest to the given point x""" |
| 78 |
if isinstance(where,FunctionSpace): |
| 79 |
self.__where=where |
| 80 |
else: |
| 81 |
self.__where=ContinuousFunction(where) |
| 82 |
self.__id=length(x-self.__where.getX()).minarg() |
| 83 |
|
| 84 |
def getValue(self,data): |
| 85 |
"""returns the value of data at the Location at a numarray object""" |
| 86 |
if isinstance(data,Data): |
| 87 |
return data |
| 88 |
else: |
| 89 |
if not data.getFunctionSpace()==self.getFunctionSpace(): |
| 90 |
raise ValueError,"function space of data obejct does not match function space of Location" |
| 91 |
else: |
| 92 |
return data.getValue(self.getId()) |
| 93 |
def getX(self): |
| 94 |
"""returns the exact coordinates of the Location""" |
| 95 |
return self.getValue(self.getFunctionSpace().getX()) |
| 96 |
|
| 97 |
def getId(self): |
| 98 |
"""returns the identifier of the location""" |
| 99 |
return self.__id |
| 100 |
|
| 101 |
def getFunctionSpace(self): |
| 102 |
"""returns the function space of the Location""" |
| 103 |
return self.__function_space |
| 104 |
|
| 105 |
def __str__(self): |
| 106 |
"""returns the coordinates of the Location as a string""" |
| 107 |
return str(self.getX()) |
| 108 |
|
| 109 |
def testProjector(domain): |
| 110 |
"""runs a few test of the Projector factory and returns the largest error plus a description of the test this error occured""" |
| 111 |
error_max=0. |
| 112 |
error_text="" |
| 113 |
x=ContinuousFunction(domain).getX() |
| 114 |
for f in [True,False]: |
| 115 |
p=Projector(domain,reduce=False,fast=f) |
| 116 |
for r in range(5): |
| 117 |
text="range %s , fast=%s"%(r,f) |
| 118 |
if r==0: |
| 119 |
td_ref=x[0] |
| 120 |
elif r==1: |
| 121 |
td_ref=x |
| 122 |
elif r==2: |
| 123 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 124 |
elif r==3: |
| 125 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 126 |
elif r==3: |
| 127 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 128 |
td=p(td_ref.interpolate(Function(domain))) |
| 129 |
er=Lsup(td-td_ref)/Lsup(td_ref) |
| 130 |
print text," error = ",er |
| 131 |
if er>error_max: |
| 132 |
error_max=er |
| 133 |
error_text=text |
| 134 |
return error_max,error_text |
| 135 |
|
| 136 |
|
| 137 |
# this should be removed later |
| 138 |
if __name__=="__main__": |
| 139 |
from esys.finley import Rectangle |
| 140 |
txt=testProjector(Rectangle(56,61)) |
| 141 |
print "test Projector: ",txt |
| 142 |
|
| 143 |
|
| 144 |
|
Parent Directory
| 
Revision Log