# Diff of /temp_trunk_copy/escript/py_src/linearPDEs.py

revision 430 by gross, Wed Jan 11 06:40:50 2006 UTC revision 596 by gross, Wed Mar 15 22:41:23 2006 UTC
# Line 448  class LinearPDE(object): Line 448  class LinearPDE(object):
448     AMG= 22     AMG= 22
449     RILU = 23     RILU = 23
450
451     __TOL=1.e-13     SMALL_TOLERANCE=1.e-13
452     __PACKAGE_KEY="package"     __PACKAGE_KEY="package"
453     __METHOD_KEY="method"     __METHOD_KEY="method"
454     __SYMMETRY_KEY="symmetric"     __SYMMETRY_KEY="symmetric"
# Line 703  class LinearPDE(object): Line 703  class LinearPDE(object):
703        else:        else:
704           A=self.getCoefficientOfGeneralPDE("A")           A=self.getCoefficientOfGeneralPDE("A")
705           if not A.isEmpty():           if not A.isEmpty():
706              tol=util.Lsup(A)*self.__TOL              tol=util.Lsup(A)*self.SMALL_TOLERANCE
707              if self.getNumSolutions()>1:              if self.getNumSolutions()>1:
708                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
709                    for j in range(self.getDim()):                    for j in range(self.getDim()):
# Line 727  class LinearPDE(object): Line 727  class LinearPDE(object):
727              if verbose: print "non-symmetric PDE because C is not present but B is"              if verbose: print "non-symmetric PDE because C is not present but B is"
728              out=False              out=False
729           elif not B.isEmpty() and not C.isEmpty():           elif not B.isEmpty() and not C.isEmpty():
730              tol=(util.Lsup(B)+util.Lsup(C))*self.__TOL/2.              tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
731              if self.getNumSolutions()>1:              if self.getNumSolutions()>1:
732                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
733                     for j in range(self.getDim()):                     for j in range(self.getDim()):
# Line 743  class LinearPDE(object): Line 743  class LinearPDE(object):
743           if self.getNumSolutions()>1:           if self.getNumSolutions()>1:
744             D=self.getCoefficientOfGeneralPDE("D")             D=self.getCoefficientOfGeneralPDE("D")
745             if not D.isEmpty():             if not D.isEmpty():
746               tol=util.Lsup(D)*self.__TOL               tol=util.Lsup(D)*self.SMALL_TOLERANCE
747               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
748                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
749                    if util.Lsup(D[i,k]-D[k,i])>tol:                    if util.Lsup(D[i,k]-D[k,i])>tol:
# Line 751  class LinearPDE(object): Line 751  class LinearPDE(object):
751                        out=False                        out=False
752             d=self.getCoefficientOfGeneralPDE("d")             d=self.getCoefficientOfGeneralPDE("d")
753             if not d.isEmpty():             if not d.isEmpty():
754               tol=util.Lsup(d)*self.__TOL               tol=util.Lsup(d)*self.SMALL_TOLERANCE
755               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
756                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
757                    if util.Lsup(d[i,k]-d[k,i])>tol:                    if util.Lsup(d[i,k]-d[k,i])>tol:
# Line 759  class LinearPDE(object): Line 759  class LinearPDE(object):
759                        out=False                        out=False
760             d_contact=self.getCoefficientOfGeneralPDE("d_contact")             d_contact=self.getCoefficientOfGeneralPDE("d_contact")
761             if not d_contact.isEmpty():             if not d_contact.isEmpty():
762               tol=util.Lsup(d_contact)*self.__TOL               tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
763               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
764                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
765                    if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:                    if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:
# Line 825  class LinearPDE(object): Line 825  class LinearPDE(object):
825         sets a new solver         sets a new solver
826
827         @param solver: sets a new solver method.         @param solver: sets a new solver method.
828         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{PRES20}, L{LUMPING}.         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{PRES20}, L{LUMPING}, L{AMG}
829         @param preconditioner: sets a new solver method.         @param preconditioner: sets a new solver method.
830         @type solver: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
831         """         """
832         if solver==None: solve=self.DEFAULT         if solver==None: solve=self.DEFAULT
833         if preconditioner==None: preconditioner=self.DEFAULT         if preconditioner==None: preconditioner=self.DEFAULT
# Line 860  class LinearPDE(object): Line 860  class LinearPDE(object):
860         elif m[0]==self.GMRES: method= "GMRES"         elif m[0]==self.GMRES: method= "GMRES"
861         elif m[0]==self.PRES20: method= "PRES20"         elif m[0]==self.PRES20: method= "PRES20"
862         elif m[0]==self.LUMPING: method= "LUMPING"         elif m[0]==self.LUMPING: method= "LUMPING"
863           elif m[0]==self.AMG: method= "AMG"
864         if m[1]==self.DEFAULT: method+="+DEFAULT"         if m[1]==self.DEFAULT: method+="+DEFAULT"
865         elif m[1]==self.JACOBI: method+= "+JACOBI"         elif m[1]==self.JACOBI: method+= "+JACOBI"
866         elif m[1]==self.ILU0: method+= "+ILU0"         elif m[1]==self.ILU0: method+= "+ILU0"
867         elif m[1]==self.ILUT: method+= "+ILUT"         elif m[1]==self.ILUT: method+= "+ILUT"
868         elif m[1]==self.SSOR: method+= "+SSOR"         elif m[1]==self.SSOR: method+= "+SSOR"
869           elif m[1]==self.AMG: method+= "+AMG"
870           elif m[1]==self.RILU: method+= "+RILU"
871         if p==self.DEFAULT: package="DEFAULT"         if p==self.DEFAULT: package="DEFAULT"
872         elif p==self.PASO: package= "PASO"         elif p==self.PASO: package= "PASO"
873         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
# Line 1529  class LinearPDE(object): Line 1532  class LinearPDE(object):
1532         if not self.__operator_is_Valid or not self.__righthandside_isValid:         if not self.__operator_is_Valid or not self.__righthandside_isValid:
1533            if self.isUsingLumping():            if self.isUsingLumping():
1534                if not self.__operator_is_Valid:                if not self.__operator_is_Valid:
1535                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1536                   if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Using coefficient A in lumped matrix can produce wrong results"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1537                   if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Using coefficient B in lumped matrix can produce wrong results"                   if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1538                   if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Using coefficient C in lumped matrix can produce wrong results"                        raise ValueError,"coefficient A in lumped matrix may not be present."
1539                   mat=self.__getNewOperator()                   if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1540                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                        raise ValueError,"coefficient A in lumped matrix may not be present."
1541                             self.getCoefficientOfGeneralPDE("A"), \                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1542                             self.getCoefficientOfGeneralPDE("B"), \                        raise ValueError,"coefficient A in lumped matrix may not be present."
1543                             self.getCoefficientOfGeneralPDE("C"), \                   D=self.getCoefficientOfGeneralPDE("D")
1544                             self.getCoefficientOfGeneralPDE("D"), \                   if not D.isEmpty():
1545                             escript.Data(), \                       if self.getNumSolutions()>1:
1546                             escript.Data(), \                          D_times_e=util.matrixmult(D,numarray.ones((self.getNumSolutions(),)))
1547                             self.getCoefficientOfGeneralPDE("d"), \                       else:
1548                             escript.Data(),\                          D_times_e=D
1549                             self.getCoefficientOfGeneralPDE("d_contact"), \                   else:
1550                             escript.Data())                      D_times_e=escript.Data()
1551                   self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))                   d=self.getCoefficientOfGeneralPDE("d")
1552                   del mat                   if not d.isEmpty():
1553                         if self.getNumSolutions()>1:
1554                            d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))
1555                         else:
1556                            d_times_e=d
1557                     else:
1558                        d_times_e=escript.Data()
1559                     d_contact=self.getCoefficientOfGeneralPDE("d_contact")
1560                     if not d_contact.isEmpty():
1561                         if self.getNumSolutions()>1:
1562                            d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),)))
1563                         else:
1564                            d_contact_times_e=d_contact
1565                     else:
1566                        d_contact_times_e=escript.Data()
1567
1568                     self.__operator=self.__getNewRightHandSide()
1570                                                  escript.Data(), \
1571                                                  D_times_e, \
1572                                                  d_times_e,\
1573                                                  d_contact_times_e)
1574                     self.__operator=1./self.__operator
1575                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
1576                   self.__operator_is_Valid=True                   self.__operator_is_Valid=True
1577                if not self.__righthandside_isValid:                if not self.__righthandside_isValid:
# Line 1811  class LameEquation(LinearPDE): Line 1836  class LameEquation(LinearPDE):
1836                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1837        self.setSymmetryOn()        self.setSymmetryOn()
1838
1839     def setValue(self,**coefficients):     def setValues(self,**coefficients):
1840       """       """
1841       sets new values to coefficients       sets new values to coefficients
1842
# Line 1834  class LameEquation(LinearPDE): Line 1859  class LameEquation(LinearPDE):
1859                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1860       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1861       """       """
1862       super(LameEquation, self).setValue(**coefficients)       super(LameEquation, self).setValues(**coefficients)
1863
1864     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1865       """       """
# Line 1948  class AdvectivePDE(LinearPDE): Line 1973  class AdvectivePDE(LinearPDE):
1973           self.__xi=xi           self.__xi=xi
1974        self.__Xi=escript.Data()        self.__Xi=escript.Data()
1975
1976     def setValue(**coefficients):     def setValue(self,**coefficients):
1977        """        """
1978        sets new values to coefficients        sets new values to coefficients
1979
# Line 2032  class AdvectivePDE(LinearPDE): Line 2057  class AdvectivePDE(LinearPDE):
2057       """       """
2058       return escript.Scalar(0.5,P.getFunctionSpace())       return escript.Scalar(0.5,P.getFunctionSpace())
2059
def __calculateXi(self,peclet_factor,flux,h):
flux=util.Lsup(flux)
if flux_max>0.:
return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)
else:
return 0.

2060     def __getXi(self):     def __getXi(self):
2061        if self.__Xi.isEmpty():        if self.__Xi.isEmpty():
2062           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2048  class AdvectivePDE(LinearPDE): Line 2066  class AdvectivePDE(LinearPDE):
2066           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))           self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2067           if not C.isEmpty() or not B.isEmpty():           if not C.isEmpty() or not B.isEmpty():
2068              if not C.isEmpty() and not B.isEmpty():              if not C.isEmpty() and not B.isEmpty():
flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2069                  if self.getNumEquations()>1:                  if self.getNumEquations()>1:
2070                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2071                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2072                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2073                           for k in range(self.getNumSolutions()):                           for k in range(self.getNumSolutions()):
2074                              for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2                              for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2
2075                          length_of_flux=util.sqrt(flux2)
2076                        # flux=C-util.reorderComponents(B,[0,2,1])                        # flux=C-util.reorderComponents(B,[0,2,1])
2077                     else:                     else:
2078                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2079                        for i in range(self.getNumEquations()):                        for i in range(self.getNumEquations()):
2080                           for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2                           for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2
2081                          length_of_flux=util.sqrt(flux2)
2082                        # flux=C-B                        # flux=C-B
2083                  else:                  else:
2084                     if self.getNumSolutions()>1:                     if self.getNumSolutions()>1:
2085                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2086                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
2087                           for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2                           for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2
2088                        # flux=C-util.reorderComponents(B,[1,0])                        # flux=C-util.reorderComponents(B,[1,0])
2089                          length_of_flux=util.sqrt(flux2)
2090                     else:                     else:
2091                        for l in range(self.getDim()): flux2+=(C[l]-B[l])**2                        length_of_flux=util.length(C-B)
#flux=C-B
length_of_flux=util.sqrt(flux2)
2092              elif C.isEmpty():              elif C.isEmpty():
2093                length_of_flux=util.length(B)                length_of_flux=util.length(B)
#flux=B
2094              else:              else:
2095                length_of_flux=util.length(C)                length_of_flux=util.length(C)
#flux=C

#length_of_flux=util.length(flux)
2096              flux_max=util.Lsup(length_of_flux)              flux_max=util.Lsup(length_of_flux)
2097              if flux_max>0.:              if flux_max>0.:
2098                 # length_of_A=util.inner(flux,util.tensormutiply(A,flux))                if A.isEmpty():
2099                 length_of_A=util.length(A)                    inv_A=1./self.SMALL_TOLERANCE
2100                 A_max=util.Lsup(length_of_A)                    peclet_number=escript.Scalar(inv_A,length_of_flux.getFunctionSpace())
2101                 if A_max>0:                    xi=self.__xi(self,peclet_number)
2102                      inv_A=1./(length_of_A+A_max*self.__TOL)                else:
2103                 else:                    # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2104                      inv_A=1./self.__TOL                    length_of_A=util.length(A)
2105                 peclet_number=length_of_flux*h/2*inv_A                    A_max=util.Lsup(length_of_A)
2106                 xi=self.__xi(peclet_number)                    if A_max>0:
2107                 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)                         inv_A=1./(length_of_A+A_max*self.SMALL_TOLERANCE)
2108                 self.trace("preclet number = %e"%util.Lsup(peclet_number))                    else:
2109                           inv_A=1./self.SMALL_TOLERANCE
2110                      peclet_number=length_of_flux*h/2*inv_A
2111                      xi=self.__xi(self,peclet_number)
2112                  self.__Xi=h*xi/(length_of_flux+flux_max*self.SMALL_TOLERANCE)
2113                  self.trace("preclet number = %e"%util.Lsup(peclet_number))
2114                else:
2115                  self.__Xi=escript.Scalar(0.,length_of_flux.getFunctionSpace())
2116        return self.__Xi        return self.__Xi
2117
2118
# Line 2115  class AdvectivePDE(LinearPDE): Line 2139  class AdvectivePDE(LinearPDE):
2139              Aout=A              Aout=A
2140           else:           else:
2141              if A.isEmpty():              if A.isEmpty():
2142                 Aout=self.createNewCoefficient("A")                 Aout=self.createCoefficientOfGeneralPDE("A")
2143              else:              else:
2144                 Aout=A[:]                 Aout=A[:]
2145              Xi=self.__getXi()              Xi=self.__getXi()
# Line 2135  class AdvectivePDE(LinearPDE): Line 2159  class AdvectivePDE(LinearPDE):
2159                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
2160                                 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)                                 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2161              else:              else:
2162                  for j in range(self.getDim()):                 if not C.isEmpty() and not B.isEmpty():
2163                     for l in range(self.getDim()):                     delta=(C-B)
2164                        if not C.isEmpty() and not B.isEmpty():                     Aout+=util.outer(Xi*delta,delta)
2165                            Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])                 elif not B.isEmpty():
2166                        elif C.isEmpty():                     Aout+=util.outer(Xi*B,B)
2167                            Aout[j,l]+=Xi*B[j]*B[l]                 elif not C.isEmpty():
2168                        else:                     Aout+=util.outer(Xi*C,C)
Aout[j,l]+=Xi*C[j]*C[l]
# if not C.isEmpty() and not B.isEmpty():
#    tmp=C-B
#    Aout=Aout+Xi*util.outer(tmp,tmp)
# elif C.isEmpty():
#    Aout=Aout+Xi*util.outer(B,B)
# else:
# Aout=Aout+Xi*util.outer(C,C)
2169           return Aout           return Aout
2170       elif name == "B" :       elif name == "B" :
2171             # return self.getCoefficient("B")
2172           B=self.getCoefficient("B")           B=self.getCoefficient("B")
2173           C=self.getCoefficient("C")           C=self.getCoefficient("C")
2174           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 2160  class AdvectivePDE(LinearPDE): Line 2177  class AdvectivePDE(LinearPDE):
2177           else:           else:
2178              Xi=self.__getXi()              Xi=self.__getXi()
2179              if B.isEmpty():              if B.isEmpty():
2180                  Bout=self.createNewCoefficient("B")                  Bout=self.createCoefficientOfGeneralPDE("B")
2181              else:              else:
2182                  Bout=B[:]                  Bout=B[:]
2183              if self.getNumEquations()>1:              if self.getNumEquations()>1:
# Line 2172  class AdvectivePDE(LinearPDE): Line 2189  class AdvectivePDE(LinearPDE):
2189                             Bout[i,j,k]+=tmp*C[p,i,j]                             Bout[i,j,k]+=tmp*C[p,i,j]
2190                             # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)                             # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2191              else:              else:
2192                 tmp=Xi*D                 Bout+=(Xi*D)*C
for j in range(self.getDim()): Bout[j]+=tmp*C[j]
# Bout=Bout+Xi*D*C
2193           return Bout           return Bout
2194       elif name == "C" :       elif name == "C" :
2195             # return self.getCoefficient("C")
2196           B=self.getCoefficient("B")           B=self.getCoefficient("B")
2197           C=self.getCoefficient("C")           C=self.getCoefficient("C")
2198           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 2185  class AdvectivePDE(LinearPDE): Line 2201  class AdvectivePDE(LinearPDE):
2201           else:           else:
2202              Xi=self.__getXi()              Xi=self.__getXi()
2203              if C.isEmpty():              if C.isEmpty():
2204                  Cout=self.createNewCoefficient("C")                  Cout=self.createCoefficientOfGeneralPDE("C")
2205              else:              else:
2206                  Cout=C[:]                  Cout=C[:]
2207              if self.getNumEquations()>1:              if self.getNumEquations()>1:
# Line 2197  class AdvectivePDE(LinearPDE): Line 2213  class AdvectivePDE(LinearPDE):
2213                                   Cout[i,k,l]+=tmp*B[p,l,i]                                   Cout[i,k,l]+=tmp*B[p,l,i]
2214                                   # Cout=Cout+Xi*B[p,l,i]*D[p,k]                                   # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2215              else:              else:
2216                 tmp=Xi*D                 Cout+=(Xi*D)*B
for j in range(self.getDim()): Cout[j]+=tmp*B[j]
# Cout=Cout+tmp*D*B
2217           return Cout           return Cout
2218       elif name == "D" :       elif name == "D" :
2219           return self.getCoefficient("D")           return self.getCoefficient("D")
2220       elif name == "X" :       elif name == "X" :
2221             # return self.getCoefficient("X")
2222           X=self.getCoefficient("X")           X=self.getCoefficient("X")
2223           Y=self.getCoefficient("Y")           Y=self.getCoefficient("Y")
2224           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 2212  class AdvectivePDE(LinearPDE): Line 2227  class AdvectivePDE(LinearPDE):
2227              Xout=X              Xout=X
2228           else:           else:
2229              if X.isEmpty():              if X.isEmpty():
2230                  Xout=self.createNewCoefficient("X")                  Xout=self.createCoefficientOfGeneralPDE("X")
2231              else:              else:
2232                  Xout=X[:]                  Xout=X[:]
2233              Xi=self.__getXi()              Xi=self.__getXi()
# Line 2231  class AdvectivePDE(LinearPDE): Line 2246  class AdvectivePDE(LinearPDE):
2246                               Xout[i,j]+=tmp*C[p,i,j]                               Xout[i,j]+=tmp*C[p,i,j]
2247                               # Xout=X_out+Xi*util.inner(Y,C,offset=1)                               # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2248              else:              else:
2249                   tmp=Xi*Y                if not C.isEmpty() and not B.isEmpty():
2250                   for j in range(self.getDim()):                  Xout+=(Xi*Y)*(C-B)
2251                      if not C.isEmpty() and not B.isEmpty():                elif C.isEmpty():
2252                         Xout[j]+=tmp*(C[j]-B[j])                  Xout-=(Xi*Y)*B
2253                         # Xout=Xout+Xi*Y*(C-B)                else:
2254                      elif C.isEmpty():                  Xout+=(Xi*Y)*C
Xout[j]-=tmp*B[j]
# Xout=Xout-Xi*Y*B
else:
Xout[j]+=tmp*C[j]
# Xout=Xout+Xi*Y*C
2255           return Xout           return Xout
2256       elif name == "Y" :       elif name == "Y" :
2257           return self.getCoefficient("Y")           return self.getCoefficient("Y")
# Line 2260  class AdvectivePDE(LinearPDE): Line 2270  class AdvectivePDE(LinearPDE):
2270       else:       else:
2271          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2272
"""
Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form

with natural boundary conditons

M{n[j]*matrixmult(k,grad(u))[j] = g- S{alpha}u }

and constraints:

M{u=r} where M{q>0}

and

M{k_bar[i,j]=k[i,j]+upwind[i]*upwind[j]}

"""

def __init__(self,domain,debug=False):
"""
initializes a new Poisson equation

@param domain: domain of the PDE
@type domain: L{Domain<escript.Domain>}
@param debug: if True debug informations are printed.

"""
self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
"k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
"v": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
"upwind": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
"alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
"g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
"r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
"q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}

def setValue(self,**coefficients):
"""
sets new values to coefficients

@param coefficients: new values assigned to coefficients
@keyword omega: value for coefficient M{S{omega}}
@type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
@keyword k: value for coefficient M{k}
@type k: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}.
@keyword v: value for coefficient M{v}
@type v: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
@keyword upwind: value for upwind term M{upwind}
@type upwind: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
@keyword f: value for right hand side M{f}
@type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
@keyword alpha: value for right hand side M{S{alpha}}
@type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
@keyword g: value for right hand side M{g}
@type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
@keyword r: prescribed values M{r} for the solution in constraints.
@type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
depending of reduced order is used for the representation of the equation.
@keyword q: mask for location of constraints
@type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
depending of reduced order is used for the representation of the equation.
@raise IllegalCoefficient: if an unknown coefficient keyword is used.
"""

def getCoefficientOfGeneralPDE(self,name):
"""
return the value of the coefficient name of the general PDE

@param name: name of the coefficient requested.
@type name: C{string}
@return: the value of the coefficient  name
@rtype: L{Data<escript.Data>}
@raise IllegalCoefficient: if name is not one of coefficients
"A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
@note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
"""
if name == "A" :
return self.getCoefficient("k")+util.outer(self.getCoefficient("upwind"),self.getCoefficient("upwind"))
elif name == "B" :
return escript.Data()
elif name == "C" :
return self.getCoefficient("v")
elif name == "D" :
return self.getCoefficient("omega")
elif name == "X" :
return escript.Data()
elif name == "Y" :
return self.getCoefficient("f")
elif name == "d" :
return self.getCoefficient("alpha")
elif name == "y" :
return self.getCoefficient("g")
elif name == "d_contact" :
return escript.Data()
elif name == "y_contact" :
return escript.Data()
elif name == "r" :
return self.getCoefficient("r")
elif name == "q" :
return self.getCoefficient("q")
else:
raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name

2273  # \$Log\$  # \$Log\$
2274  # Revision 1.14  2005/09/22 01:54:57  jgs  # Revision 1.14  2005/09/22 01:54:57  jgs
2275  # Merge of development branch dev-02 back to main trunk on 2005-09-22  # Merge of development branch dev-02 back to main trunk on 2005-09-22

Legend:
 Removed from v.430 changed lines Added in v.596