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