# Diff of /trunk/finley/test/python/tp.py

revision 1368 by gross, Wed Dec 19 03:45:20 2007 UTC revision 1370 by gross, Wed Jan 2 09:21:43 2008 UTC
# Line 1  Line 1
1  from esys.escript import *  from esys.escript import *
2
3  class TransportPDE(object):  class TransportPDE(object):
4      def __init__(domain,num_equations=1,theta=0.,dt_max=-1.,trace=True):       def __init__(self,domain,num_equations=1,theta=0.,dt_max=-1.,trace=True):
5          self.__domain=domain          self.__domain=domain
6          self.__num_equations=num_equations          self.__num_equations=num_equations
7          self.__theta=theta          self.__theta=theta
8          self.__dt_max=dt_max          self.__dt_max=dt_max
9          self.__transport_problem=None          self.__transport_problem=None
10        self.__trace=trace
11        self.__matrix_type=0
12         def trace(self,text):
13             if self.__trace: print text
14
15       def getDomain(self):       def getDomain(self):
16          return self.__domain          return self.__domain
# Line 15  class TransportPDE(object): Line 19  class TransportPDE(object):
19       def getDt_max(self):       def getDt_max(self):
20          return self.__dt_max          return self.__dt_max
21       def getNumEquations(self):       def getNumEquations(self):
22          return self.__getNumEquations          return self.__num_equations
23         def reduced(self):
24             return False
25       def getFunctionSpace(self):       def getFunctionSpace(self):
26          if self.redueced():          if self.reduced():
27             return ReducedSolution(self.getDomain())             return ReducedSolution(self.getDomain())
28          else:          else:
29             return Solution(self.getDomain())             return Solution(self.getDomain())
# Line 35  class TransportPDE(object): Line 41  class TransportPDE(object):
41                                 self.getFunctionSpace(), \                                 self.getFunctionSpace(), \
42                                 self.__matrix_type)                                 self.__matrix_type)
43
44       def setValue(M=Data(),A=Data(),B=Data(),C=Data(),D=Data(),X=Data(),Y=Data()):       def setValue(self,M=Data(),A=Data(),B=Data(),C=Data(),D=Data(),X=Data(),Y=Data(),
45                  d=Data(),y=Data(),d_contact=Data(),y_contact=Data()):
46           self.__transport_problem=self.__getNewTransportProblem()           self.__transport_problem=self.__getNewTransportProblem()
47           self.getDomain().       if self.getNumEquations() ==1 :
48            self.__source=Data(0.0,(),self.getFunctionSpace())
49       def solve(X=Data(),Y=Data()):           else:
50                 self.__source=Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())
52                     self.__transport_problem,
53                 self.__source,
54                 M,A,B,C,D,X,Y,d,y,d_contact,y_contact)
55
56         def setInitialSolution(self,u):
57             self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace()))
58         def solve(self,dt):
59           return self.__transport_problem.solve(self.__source,dt,{})
60    from esys.finley import Rectangle
61
62    dom=Rectangle(10,10)
63    fc=TransportPDE(dom,num_equations=1,theta=0.0,dt_max=0.)
64    fc.setValue(M=Scalar(1.,Function(dom)),C=Scalar(1.,Function(dom))*[-1.,0.])
65    x=dom.getX()
66    u=whereNonPositive(abs(x[0]-0.3)-0.1)*whereNonPositive(abs(x[1]-0.3)-0.1)
67    c=0
68    saveVTK("u.%s.xml"%c,u=u)
69    fc.setInitialSolution(u)
70    dt=0.301
71    t=0.
72    while t<dt:
73        print "time step t=",t+dt
74        u=fc.solve(dt)
75        print "range u",inf(u),sup(u)
76        c+=1
77        saveVTK("u.%s.xml"%c,u=u)
78        t+=dt

Legend:
 Removed from v.1368 changed lines Added in v.1370