--- trunk/escript/py_src/linearPDEs.py 2006/01/10 04:10:39 425 +++ trunk/escript/py_src/linearPDEs.py 2007/04/24 08:55:04 1118 @@ -1,14 +1,4 @@ # $Id$ - -# -# COPYRIGHT ACcESS 2004 - All Rights Reserved -# -# This software is the property of ACcESS. No part of this code -# may be copied in any form or by any means without the expressed written -# consent of ACcESS. Copying, use or modification of this software -# by any unauthorised person is illegal unless that -# person has a software license agreement with ACcESS. -# """ The module provides an interface to define and solve linear partial differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any @@ -17,12 +7,13 @@ The general interface is provided through the L{LinearPDE} class. The L{AdvectivePDE} which is derived from the L{LinearPDE} class provides an interface to PDE dominated by its advective terms. The L{Poisson}, -L{Helmholtz}, L{LameEquation}, L{AdvectionDiffusion} +L{Helmholtz}, L{LameEquation}, L{AdvectivePDE} classs which are also derived form the L{LinearPDE} class should be used to define of solve these sepecial PDEs. @var __author__: name of author -@var __licence__: licence agreement +@var __copyright__: copyrights +@var __license__: licence agreement @var __url__: url entry point on documentation @var __version__: version @var __date__: date of the version @@ -33,8 +24,12 @@ import numarray __author__="Lutz Gross, l.gross@uq.edu.au" -__licence__="contact: esys@access.uq.edu.au" -__url__="http://www.iservo.edu.au/esys/escript" +__copyright__=""" Copyright (c) 2006 by ACcESS MNRF + http://www.access.edu.au + Primary Business: Queensland, Australia""" +__license__="""Licensed under the Open Software License version 3.0 + http://www.opensource.org/licenses/osl-3.0.php""" +__url__="http://www.iservo.edu.au/esys" __version__="$Revision$" __date__="$Date$" @@ -43,16 +38,24 @@ """ raised if an illegal coefficient of the general ar particular PDE is requested. """ + pass class IllegalCoefficientValue(ValueError): """ raised if an incorrect value for a coefficient is used. """ + pass + +class IllegalCoefficientFunctionSpace(ValueError): + """ + raised if an incorrect function space for a coefficient is used. + """ class UndefinedPDEError(ValueError): """ raised if a PDE is not fully defined yet. """ + pass class PDECoefficient(object): """ @@ -61,6 +64,9 @@ @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain + @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the interior of the PDE domain using a reduced integration order + @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the boundary of the PDE domain using a reduced integration order + @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact region within the PDE domain using a reduced integration order @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations @@ -82,13 +88,17 @@ OPERATOR=10 RIGHTHANDSIDE=11 BOTH=12 + INTERIOR_REDUCED=13 + BOUNDARY_REDUCED=14 + CONTACT_REDUCED=15 - def __init__(self,where,pattern,altering): + def __init__(self, where, pattern, altering): """ Initialise a PDE Coefficient type @param where: describes where the coefficient lives - @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED} + @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED}, + L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED}, L{CONTACT_REDUCED}. @param pattern: describes the shape of the coefficient and how the shape is build for a given spatial dimension and numbers of equation and solution in then PDE. For instance, (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which @@ -99,7 +109,8 @@ @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM} @param altering: indicates what part of the PDE is altered if the coefficiennt is altered @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH} - + @param reduced: indicates if reduced + @type reduced: C{bool} """ super(PDECoefficient, self).__init__() self.what=where @@ -120,18 +131,24 @@ @param domain: domain on which the PDE uses the coefficient @type domain: L{Domain} @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation - @type domain: C{bool} + @type reducedEquationOrder: C{bool} @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution - @type domain: C{bool} + @type reducedSolutionOrder: C{bool} @return: L{FunctionSpace} of the coefficient @rtype: L{FunctionSpace} """ if self.what==self.INTERIOR: return escript.Function(domain) + elif self.what==self.INTERIOR_REDUCED: + return escript.ReducedFunction(domain) elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain) + elif self.what==self.BOUNDARY_REDUCED: + return escript.ReducedFunctionOnBoundary(domain) elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain) + elif self.what==self.CONTACT_REDUCED: + return escript.ReducedFunctionOnContactZero(domain) elif self.what==self.SOLUTION: if reducedEquationOrder and reducedSolutionOrder: return escript.ReducedSolution(domain) @@ -160,21 +177,23 @@ @param numSolutions: number of components of the PDE solution @type numSolutions: C{int} @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation - @type domain: C{bool} + @type reducedEquationOrder: C{bool} @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution - @type domain: C{bool} + @type reducedSolutionOrder: C{bool} @param newValue: number of components of the PDE solution @type newValue: any object that can be converted into a L{Data} object with the appropriate shape and L{FunctionSpace} @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient + @raise IllegalCoefficientFunctionSpace: if unable to interploate value to appropriate function space """ if newValue==None: newValue=escript.Data() elif isinstance(newValue,escript.Data): if not newValue.isEmpty(): - try: - newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder)) - except: - raise IllegalCoefficientValue,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain) + if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder): + try: + newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder)) + except: + raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain) else: newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder)) if not newValue.isEmpty(): @@ -318,21 +337,22 @@ For a single PDE with a solution with a single component the linear PDE is defined in the following form: - M{-grad(A[j,l]*grad(u)[l]+B[j]u)[j]+C[l]*grad(u)[l]+D*u =-grad(X)[j,j]+Y} + M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)} + where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention, ie. summation over indexes appearing twice in a term of a sum is performed, is used. The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data} objects in the - L{Function} on the PDE or objects that can be converted into such L{Data} objects. - M{A} is a rank two, M{B}, M{C} and M{X} are rank one and M{D} and M{Y} are scalar. + L{Function} and the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data} objects in the + L{ReducedFunction}. It is also allowd to use objects that can be converted into + such L{Data} objects. M{A} and M{A_reduced} are rank two, M{B_reduced}, M{C_reduced}, M{X_reduced} + M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced} and M{Y_reduced} are scalar. The following natural boundary conditions are considered: - M{n[j]*(A[i,j]*grad(u)[l]+B[j]*u)+d*u=n[j]*X[j]+y} + M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y} - where M{n} is the outer normal field calculated by L{getNormal} of L{FunctionOnBoundary}. - Notice that the coefficients M{A}, M{B} and M{X} are defined in the PDE. The coefficients M{d} and M{y} are - each a scalar in the L{FunctionOnBoundary}. + where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the PDE. The coefficients M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary} and the coefficients M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary}. Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form @@ -344,20 +364,21 @@ The PDE is symmetrical if - M{A[i,j]=A[j,i]} and M{B[j]=C[j]} + M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and M{B_reduced[j]=C_reduced[j] For a system of PDEs and a solution with several components the PDE has the form - M{-grad(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])[j]+C[i,k,l]*grad(u[k])[l]+D[i,k]*u[k] =-grad(X[i,j])[j]+Y[i] } + M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] } - M{A} is a ramk four, M{B} and M{C} are each a rank three, M{D} and M{X} are each a rank two and M{Y} is a rank one. + M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one. The natural boundary conditions take the form: - M{n[j]*(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])+d[i,k]*u[k]=n[j]*X[i,j]+y[i]} + M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]} - The coefficient M{d} is a rank two and M{y} is a rank one both in the L{FunctionOnBoundary}. Constraints take the form + The coefficient M{d} is a rank two and M{y} is a rank one both in the L{FunctionOnBoundary}. Constraints take the form and the coefficients M{d_reduced} is a rank two and M{y_reduced} is a rank one both in the L{ReducedFunctionOnBoundary}. + Constraints take the form M{u[i]=r[i]} where M{q[i]>0} @@ -366,36 +387,40 @@ The system of PDEs is symmetrical if - M{A[i,j,k,l]=A[k,l,i,j]} + - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]} - M{B[i,j,k]=C[k,i,j]} + - M{B_reduced[i,j,k]=C_reduced[k,i,j]} - M{D[i,k]=D[i,k]} + - M{D_reduced[i,k]=D_reduced[i,k]} - M{d[i,k]=d[k,i]} + - M{d_reduced[i,k]=d_reduced[k,i]} L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution defined as - M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k]-X[i,j]} + M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]} For the case of single solution component and single PDE M{J} is defined - M{J_{j}=A[i,j]*grad(u)[j]+B[i]*u-X[i]} + M{J_{j}=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]} In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1 calculated from L{getNormal} of L{FunctionOnContactZero}. For a system of PDEs the contact condition takes the form - M{n[j]*J0[i,j]=n[j]*J1[i,j]=y_contact[i]- d_contact[i,k]*jump(u)[k]} + M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]} where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by L{jump}. The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero} or L{FunctionOnContactOne}. + The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero} or L{ReducedFunctionOnContactOne}. In case of a single PDE and a single component solution the contact condition takes the form - M{n[j]*J0_{j}=n[j]*J1_{j}=y_contact-d_contact*jump(u)} + M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)} - In this case the the coefficient M{d_contact} and M{y_contact} are eaach scalar - both in the L{FunctionOnContactZero} or L{FunctionOnContactOne}. + In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero} or L{FunctionOnContactOne} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero} or L{ReducedFunctionOnContactOne} @cvar DEFAULT: The default method used to solve the system of linear equations @cvar DIRECT: The direct solver based on LDU factorization @@ -419,6 +444,8 @@ @cvar MKL: Intel's MKL solver library @cvar UMFPACK: the UMFPACK library @cvar ITERATIVE: The default iterative solver + @cvar AMG: algebraic multi grid + @cvar RILU: recursive ILU """ DEFAULT= 0 @@ -443,8 +470,10 @@ UMFPACK= 16 ITERATIVE= 20 PASO= 21 + AMG= 22 + RILU = 23 - __TOL=1.e-13 + SMALL_TOLERANCE=1.e-13 __PACKAGE_KEY="package" __METHOD_KEY="method" __SYMMETRY_KEY="symmetric" @@ -480,6 +509,16 @@ "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR), "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), + "A_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR), + "B_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR), + "C_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR), + "D_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR), + "X_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE), + "Y_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), + "d_reduced" : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR), + "y_reduced" : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), + "d_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR), + "y_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.RIGHTHANDSIDE), "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.BOTH)} @@ -665,9 +704,9 @@ @rtype: L{Data} """ if u==None: - return self.getOperator()*self.getSolution() + return self.getOperator()*self.getSolution() else: - self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution()) + return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution()) def getResidual(self,u=None): """ @@ -699,7 +738,7 @@ else: A=self.getCoefficientOfGeneralPDE("A") if not A.isEmpty(): - tol=util.Lsup(A)*self.__TOL + tol=util.Lsup(A)*self.SMALL_TOLERANCE if self.getNumSolutions()>1: for i in range(self.getNumEquations()): for j in range(self.getDim()): @@ -723,7 +762,7 @@ if verbose: print "non-symmetric PDE because C is not present but B is" out=False elif not B.isEmpty() and not C.isEmpty(): - tol=(util.Lsup(B)+util.Lsup(C))*self.__TOL/2. + tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2. if self.getNumSolutions()>1: for i in range(self.getNumEquations()): for j in range(self.getDim()): @@ -739,7 +778,7 @@ if self.getNumSolutions()>1: D=self.getCoefficientOfGeneralPDE("D") if not D.isEmpty(): - tol=util.Lsup(D)*self.__TOL + tol=util.Lsup(D)*self.SMALL_TOLERANCE for i in range(self.getNumEquations()): for k in range(self.getNumSolutions()): if util.Lsup(D[i,k]-D[k,i])>tol: @@ -747,7 +786,7 @@ out=False d=self.getCoefficientOfGeneralPDE("d") if not d.isEmpty(): - tol=util.Lsup(d)*self.__TOL + tol=util.Lsup(d)*self.SMALL_TOLERANCE for i in range(self.getNumEquations()): for k in range(self.getNumSolutions()): if util.Lsup(d[i,k]-d[k,i])>tol: @@ -755,12 +794,77 @@ out=False d_contact=self.getCoefficientOfGeneralPDE("d_contact") if not d_contact.isEmpty(): - tol=util.Lsup(d_contact)*self.__TOL + tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE for i in range(self.getNumEquations()): for k in range(self.getNumSolutions()): if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol: if verbose: print "non-symmetric PDE because d_contact[%d,%d]!=d_contact[%d,%d]"%(i,k,k,i) out=False + # and now the reduced coefficients + A_reduced=self.getCoefficientOfGeneralPDE("A_reduced") + if not A_reduced.isEmpty(): + tol=util.Lsup(A_reduced)*self.SMALL_TOLERANCE + if self.getNumSolutions()>1: + for i in range(self.getNumEquations()): + for j in range(self.getDim()): + for k in range(self.getNumSolutions()): + for l in range(self.getDim()): + if util.Lsup(A_reduced[i,j,k,l]-A_reduced[k,l,i,j])>tol: + if verbose: print "non-symmetric PDE because A_reduced[%d,%d,%d,%d]!=A_reduced[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j) + out=False + else: + for j in range(self.getDim()): + for l in range(self.getDim()): + if util.Lsup(A_reduced[j,l]-A_reduced[l,j])>tol: + if verbose: print "non-symmetric PDE because A_reduced[%d,%d]!=A_reduced[%d,%d]"%(j,l,l,j) + out=False + B_reduced=self.getCoefficientOfGeneralPDE("B_reduced") + C_reduced=self.getCoefficientOfGeneralPDE("C_reduced") + if B_reduced.isEmpty() and not C_reduced.isEmpty(): + if verbose: print "non-symmetric PDE because B_reduced is not present but C_reduced is" + out=False + elif not B_reduced.isEmpty() and C_reduced.isEmpty(): + if verbose: print "non-symmetric PDE because C_reduced is not present but B_reduced is" + out=False + elif not B_reduced.isEmpty() and not C_reduced.isEmpty(): + tol=(util.Lsup(B_reduced)+util.Lsup(C_reduced))*self.SMALL_TOLERANCE/2. + if self.getNumSolutions()>1: + for i in range(self.getNumEquations()): + for j in range(self.getDim()): + for k in range(self.getNumSolutions()): + if util.Lsup(B_reduced[i,j,k]-C_reduced[k,i,j])>tol: + if verbose: print "non-symmetric PDE because B_reduced[%d,%d,%d]!=C_reduced[%d,%d,%d]"%(i,j,k,k,i,j) + out=False + else: + for j in range(self.getDim()): + if util.Lsup(B_reduced[j]-C_reduced[j])>tol: + if verbose: print "non-symmetric PDE because B_reduced[%d]!=C_reduced[%d]"%(j,j) + out=False + if self.getNumSolutions()>1: + D_reduced=self.getCoefficientOfGeneralPDE("D_reduced") + if not D_reduced.isEmpty(): + tol=util.Lsup(D_reduced)*self.SMALL_TOLERANCE + for i in range(self.getNumEquations()): + for k in range(self.getNumSolutions()): + if util.Lsup(D_reduced[i,k]-D_reduced[k,i])>tol: + if verbose: print "non-symmetric PDE because D_reduced[%d,%d]!=D_reduced[%d,%d]"%(i,k,k,i) + out=False + d_reduced=self.getCoefficientOfGeneralPDE("d_reduced") + if not d_reduced.isEmpty(): + tol=util.Lsup(d_reduced)*self.SMALL_TOLERANCE + for i in range(self.getNumEquations()): + for k in range(self.getNumSolutions()): + if util.Lsup(d_reduced[i,k]-d_reduced[k,i])>tol: + if verbose: print "non-symmetric PDE because d_reduced[%d,%d]!=d_reduced[%d,%d]"%(i,k,k,i) + out=False + d_contact_reduced=self.getCoefficientOfGeneralPDE("d_contact_reduced") + if not d_contact_reduced.isEmpty(): + tol=util.Lsup(d_contact_reduced)*self.SMALL_TOLERANCE + for i in range(self.getNumEquations()): + for k in range(self.getNumSolutions()): + if util.Lsup(d_contact_reduced[i,k]-d_contact_reduced[k,i])>tol: + if verbose: print "non-symmetric PDE because d_contact_reduced[%d,%d]!=d_contact_reduced[%d,%d]"%(i,k,k,i) + out=False return out def getSolution(self,**options): @@ -800,11 +904,11 @@ """ returns the flux M{J} for a given M{u} - M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]u[k]-X[i,j]} + M{J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]} or - M{J[j]=A[i,j]*grad(u)[l]+B[j]u-X[j]} + M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]} @param u: argument in the flux. If u is not present or equals L{None} the current solution is used. @type u: L{Data} or None @@ -812,7 +916,12 @@ @rtype: L{Data} """ if u==None: u=self.getSolution() - return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u))+util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u)-util.self.getCoefficientOfGeneralPDE("X") + return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u,Funtion(self.getDomain))) \ + +util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u) \ + -util.self.getCoefficientOfGeneralPDE("X") \ + +util.tensormult(self.getCoefficientOfGeneralPDE("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \ + +util.matrixmult(self.getCoefficientOfGeneralPDE("B_reduced"),u) \ + -util.self.getCoefficientOfGeneralPDE("X_reduced") # ============================================================================= # solver settings: # ============================================================================= @@ -821,9 +930,9 @@ sets a new solver @param solver: sets a new solver method. - @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} @param preconditioner: sets a new solver method. - @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} """ if solver==None: solve=self.DEFAULT if preconditioner==None: preconditioner=self.DEFAULT @@ -856,11 +965,14 @@ elif m[0]==self.GMRES: method= "GMRES" elif m[0]==self.PRES20: method= "PRES20" elif m[0]==self.LUMPING: method= "LUMPING" + elif m[0]==self.AMG: method= "AMG" if m[1]==self.DEFAULT: method+="+DEFAULT" elif m[1]==self.JACOBI: method+= "+JACOBI" elif m[1]==self.ILU0: method+= "+ILU0" elif m[1]==self.ILUT: method+= "+ILUT" elif m[1]==self.SSOR: method+= "+SSOR" + elif m[1]==self.AMG: method+= "+AMG" + elif m[1]==self.RILU: method+= "+RILU" if p==self.DEFAULT: package="DEFAULT" elif p==self.PASO: package= "PASO" elif p==self.MKL: package= "MKL" @@ -883,12 +995,12 @@ """ sets a new solver package - @param solver: sets a new solver method. - @type solver: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMLPACK} + @param package: sets a new solver method. + @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK} """ if package==None: package=self.DEFAULT if not package==self.getSolverPackage(): - self.__solver_method=solver + self.__solver_package=package self.__checkMatrixType() self.trace("New solver is %s"%self.getSolverMethodName()) @@ -921,10 +1033,10 @@ @param tol: new tolerance for the solver. If the tol is lower then the current tolerence the system will be resolved. @type tol: positive C{float} - @raise ValueException: if tolerance is not positive. + @raise ValueError: if tolerance is not positive. """ if not tol>0: - raise ValueException,"Tolerance as to be positive" + raise ValueError,"Tolerance as to be positive" if tol} @raise IllegalCoefficient: if name is not one of coefficients - M{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}. + M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, + M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}. """ if self.hasCoefficientOfGeneralPDE(name): return self.getCoefficient(name) @@ -1283,7 +1396,8 @@ @return: a coefficient name initialized to 0. @rtype: L{Data} @raise IllegalCoefficient: if name is not one of coefficients - M{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}. + M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, + M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}. """ if self.hasCoefficientOfGeneralPDE(name): return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name)) @@ -1299,7 +1413,8 @@ @return: the function space to be used for coefficient name @rtype: L{FunctionSpace} @raise IllegalCoefficient: if name is not one of coefficients - M{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}. + M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, + M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}. """ if self.hasCoefficientOfGeneralPDE(name): return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain()) @@ -1315,7 +1430,8 @@ @return: the shape of the coefficient name @rtype: C{tuple} of C{int} @raise IllegalCoefficient: if name is not one of coefficients - M{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}. + M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, + M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}. """ if self.hasCoefficientOfGeneralPDE(name): return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions()) @@ -1445,26 +1561,42 @@ @param coefficients: new values assigned to coefficients @keyword A: value for coefficient A. @type A: any type that can be casted to L{Data} object on L{Function}. + @keyword A_reduced: value for coefficient A_reduced. + @type A_reduced: any type that can be casted to L{Data} object on L{ReducedFunction}. @keyword B: value for coefficient B @type B: any type that can be casted to L{Data} object on L{Function}. + @keyword B_reduced: value for coefficient B_reduced + @type B_reduced: any type that can be casted to L{Data} object on L{ReducedFunction}. @keyword C: value for coefficient C @type C: any type that can be casted to L{Data} object on L{Function}. + @keyword C_reduced: value for coefficient C_reduced + @type C_reduced: any type that can be casted to L{Data} object on L{ReducedFunction}. @keyword D: value for coefficient D @type D: any type that can be casted to L{Data} object on L{Function}. + @keyword D_reduced: value for coefficient D_reduced + @type D_reduced: any type that can be casted to L{Data} object on L{ReducedFunction}. @keyword X: value for coefficient X @type X: any type that can be casted to L{Data} object on L{Function}. + @keyword X_reduced: value for coefficient X_reduced + @type X_reduced: any type that can be casted to L{Data} object on L{ReducedFunction}. @keyword Y: value for coefficient Y @type Y: any type that can be casted to L{Data} object on L{Function}. + @keyword Y_reduced: value for coefficient Y_reduced + @type Y_reduced: any type that can be casted to L{Data} object on L{ReducedFunction}. @keyword d: value for coefficient d @type d: any type that can be casted to L{Data} object on L{FunctionOnBoundary}. + @keyword d_reduced: value for coefficient d_reduced + @type d_reduced: any type that can be casted to L{Data} object on L{ReducedFunctionOnBoundary}. @keyword y: value for coefficient y @type y: any type that can be casted to L{Data} object on L{FunctionOnBoundary}. @keyword d_contact: value for coefficient d_contact - @type d_contact: any type that can be casted to L{Data} object on L{FunctionOnContactOne}. - or L{FunctionOnContactZero}. + @type d_contact: any type that can be casted to L{Data} object on L{FunctionOnContactOne} or L{FunctionOnContactZero}. + @keyword d_contact_reduced: value for coefficient d_contact_reduced + @type d_contact_reduced: any type that can be casted to L{Data} object on L{ReducedFunctionOnContactOne} or L{ReducedFunctionOnContactZero}. @keyword y_contact: value for coefficient y_contact - @type y_contact: any type that can be casted to L{Data} object on L{FunctionOnContactOne}. - or L{FunctionOnContactZero}. + @type y_contact: any type that can be casted to L{Data} object on L{FunctionOnContactOne} or L{FunctionOnContactZero}. + @keyword y_contact_reduced: value for coefficient y_contact_reduced + @type y_contact_reduced: any type that can be casted to L{Data} object on L{ReducedFunctionOnContactOne} or L{ReducedFunctionOnContactZero}. @keyword r: values prescribed to the solution at the locations of constraints @type r: any type that can be casted to L{Data} object on L{Solution} or L{ReducedSolution} depending of reduced order is used for the solution. @@ -1499,11 +1631,27 @@ # now we check the shape of the coefficient if numEquations and numSolutions are set: for i,d in coefficients.iteritems(): try: - self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),self.reduceEquationOrder(),self.reduceSolutionOrder(),d) + self.COEFFICIENTS[i].setValue(self.getDomain(), + self.getNumEquations(),self.getNumSolutions(), + self.reduceEquationOrder(),self.reduceSolutionOrder(),d) + self.alteredCoefficient(i) + except IllegalCoefficientFunctionSpace,m: + # if the function space is wrong then we try the reduced version: + i_red=i+"_reduced" + if (not i_red in coefficients.keys()) and i_red in self.COEFFICIENTS.keys(): + try: + self.COEFFICIENTS[i_red].setValue(self.getDomain(), + self.getNumEquations(),self.getNumSolutions(), + self.reduceEquationOrder(),self.reduceSolutionOrder(),d) + self.alteredCoefficient(i_red) + except IllegalCoefficientValue,m: + raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m)) + except IllegalCoefficientFunctionSpace,m: + raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m)) + else: + raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m)) except IllegalCoefficientValue,m: raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m)) - self.alteredCoefficient(i) - self.__altered_coefficients=True # check if the systrem is inhomogeneous: if len(coefficients)>0 and not self.isUsingLumping(): @@ -1511,7 +1659,7 @@ r=self.getCoefficientOfGeneralPDE("r") homogeneous_constraint=True if not q.isEmpty() and not r.isEmpty(): - if util.Lsup(q*r)>=1.e-13*util.Lsup(r): + if util.Lsup(q*r)>0.: self.trace("Inhomogeneous constraint detected.") self.__invalidateSystem() @@ -1525,24 +1673,88 @@ if not self.__operator_is_Valid or not self.__righthandside_isValid: if self.isUsingLumping(): if not self.__operator_is_Valid: - if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns" - if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Using coefficient A in lumped matrix can produce wrong results" - if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Using coefficient B in lumped matrix can produce wrong results" - if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Using coefficient C in lumped matrix can produce wrong results" - mat=self.__getNewOperator() - self.getDomain().addPDEToSystem(mat,escript.Data(), \ - self.getCoefficientOfGeneralPDE("A"), \ - self.getCoefficientOfGeneralPDE("B"), \ - self.getCoefficientOfGeneralPDE("C"), \ - self.getCoefficientOfGeneralPDE("D"), \ - escript.Data(), \ - escript.Data(), \ - self.getCoefficientOfGeneralPDE("d"), \ - escript.Data(),\ - self.getCoefficientOfGeneralPDE("d_contact"), \ - escript.Data()) - self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)) - del mat + if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): + raise TypeError,"Lumped matrix requires same order for equations and unknowns" + if not self.getCoefficientOfGeneralPDE("A").isEmpty(): + raise ValueError,"coefficient A in lumped matrix may not be present." + if not self.getCoefficientOfGeneralPDE("B").isEmpty(): + raise ValueError,"coefficient B in lumped matrix may not be present." + if not self.getCoefficientOfGeneralPDE("C").isEmpty(): + raise ValueError,"coefficient C in lumped matrix may not be present." + if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty(): + raise ValueError,"coefficient A_reduced in lumped matrix may not be present." + if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty(): + raise ValueError,"coefficient B_reduced in lumped matrix may not be present." + if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty(): + raise ValueError,"coefficient C_reduced in lumped matrix may not be present." + D=self.getCoefficientOfGeneralPDE("D") + if not D.isEmpty(): + if self.getNumSolutions()>1: + D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),))) + else: + D_times_e=D + else: + D_times_e=escript.Data() + d=self.getCoefficientOfGeneralPDE("d") + if not d.isEmpty(): + if self.getNumSolutions()>1: + d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),))) + else: + d_times_e=d + else: + d_times_e=escript.Data() + d_contact=self.getCoefficientOfGeneralPDE("d_contact") + if not d_contact.isEmpty(): + if self.getNumSolutions()>1: + d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),))) + else: + d_contact_times_e=d_contact + else: + d_contact_times_e=escript.Data() + + self.__operator=self.__getNewRightHandSide() + self.getDomain().addPDEToRHS(self.__operator, \ + escript.Data(), \ + D_times_e, \ + d_times_e,\ + d_contact_times_e) + D_reduced=self.getCoefficientOfGeneralPDE("D_reduced") + if not D_reduced.isEmpty(): + if self.getNumSolutions()>1: + D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),))) + else: + D_reduced_times_e=D_reduced + else: + D_reduced_times_e=escript.Data() + d_reduced=self.getCoefficientOfGeneralPDE("d_reduced") + if not d_reduced.isEmpty(): + if self.getNumSolutions()>1: + d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),))) + else: + d_reduced_times_e=d_reduced + else: + d_reduced_times_e=escript.Data() + d_contact_reduced=self.getCoefficientOfGeneralPDE("d_contact_reduced") + if not d_contact_reduced.isEmpty(): + if self.getNumSolutions()>1: + d_contact_reduced_times_e=util.matrixmult(d_contact_reduced,numarray.ones((self.getNumSolutions(),))) + else: + d_contact_reduced_times_e=d_contact_reduced + else: + d_contact_reduced_times_e=escript.Data() + + self.__operator=self.__getNewRightHandSide() + self.getDomain().addPDEToRHS(self.__operator, \ + escript.Data(), \ + D_times_e, \ + d_times_e,\ + d_contact_times_e) + self.getDomain().addPDEToRHS(self.__operator, \ + escript.Data(), \ + D_reduced_times_e, \ + d_reduced_times_e,\ + d_contact_reduced_times_e) + self.__operator=1./self.__operator self.trace("New lumped operator has been built.") self.__operator_is_Valid=True if not self.__righthandside_isValid: @@ -1551,6 +1763,11 @@ self.getCoefficientOfGeneralPDE("Y"),\ self.getCoefficientOfGeneralPDE("y"),\ self.getCoefficientOfGeneralPDE("y_contact")) + self.getDomain().addPDEToRHS(self.__righthandside, \ + self.getCoefficientOfGeneralPDE("X_reduced"), \ + self.getCoefficientOfGeneralPDE("Y_reduced"),\ + self.getCoefficientOfGeneralPDE("y_reduced"),\ + self.getCoefficientOfGeneralPDE("y_contact_reduced")) self.trace("New right hand side as been built.") self.__righthandside_isValid=True else: @@ -1566,6 +1783,17 @@ self.getCoefficientOfGeneralPDE("y"), \ self.getCoefficientOfGeneralPDE("d_contact"), \ self.getCoefficientOfGeneralPDE("y_contact")) + self.getDomain().addPDEToSystem(self.__operator,self.__righthandside, \ + self.getCoefficientOfGeneralPDE("A_reduced"), \ + self.getCoefficientOfGeneralPDE("B_reduced"), \ + self.getCoefficientOfGeneralPDE("C_reduced"), \ + self.getCoefficientOfGeneralPDE("D_reduced"), \ + self.getCoefficientOfGeneralPDE("X_reduced"), \ + self.getCoefficientOfGeneralPDE("Y_reduced"), \ + self.getCoefficientOfGeneralPDE("d_reduced"), \ + self.getCoefficientOfGeneralPDE("y_reduced"), \ + self.getCoefficientOfGeneralPDE("d_contact_reduced"), \ + self.getCoefficientOfGeneralPDE("y_contact_reduced")) self.__applyConstraint() self.__righthandside=self.copyConstraint(self.__righthandside) self.trace("New system has been built.") @@ -1577,6 +1805,11 @@ self.getCoefficientOfGeneralPDE("Y"),\ self.getCoefficientOfGeneralPDE("y"),\ self.getCoefficientOfGeneralPDE("y_contact")) + self.getDomain().addPDEToRHS(self.__righthandside, \ + self.getCoefficientOfGeneralPDE("X_reduced"), \ + self.getCoefficientOfGeneralPDE("Y_reduced"),\ + self.getCoefficientOfGeneralPDE("y_reduced"),\ + self.getCoefficientOfGeneralPDE("y_contact_reduced")) self.__righthandside=self.copyConstraint(self.__righthandside) self.trace("New right hand side has been built.") self.__righthandside_isValid=True @@ -1592,6 +1825,17 @@ escript.Data(),\ self.getCoefficientOfGeneralPDE("d_contact"), \ escript.Data()) + self.getDomain().addPDEToSystem(self.__operator,escript.Data(), \ + self.getCoefficientOfGeneralPDE("A_reduced"), \ + self.getCoefficientOfGeneralPDE("B_reduced"), \ + self.getCoefficientOfGeneralPDE("C_reduced"), \ + self.getCoefficientOfGeneralPDE("D_reduced"), \ + escript.Data(), \ + escript.Data(), \ + self.getCoefficientOfGeneralPDE("d_reduced"), \ + escript.Data(),\ + self.getCoefficientOfGeneralPDE("d_contact_reduced"), \ + escript.Data()) self.__applyConstraint() self.trace("New operator has been built.") self.__operator_is_Valid=True @@ -1625,7 +1869,8 @@ """ super(Poisson, self).__init__(domain,1,1,debug) self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), - "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)} + "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), + "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)} self.setSymmetryOn() def setValue(self,**coefficients): @@ -1673,6 +1918,26 @@ return escript.Data() elif name == "y_contact" : return escript.Data() + elif name == "A_reduced" : + return escript.Data() + elif name == "B_reduced" : + return escript.Data() + elif name == "C_reduced" : + return escript.Data() + elif name == "D_reduced" : + return escript.Data() + elif name == "X_reduced" : + return escript.Data() + elif name == "Y_reduced" : + return self.getCoefficient("f_reduced") + elif name == "d_reduced" : + return escript.Data() + elif name == "y_reduced" : + return escript.Data() + elif name == "d_contact_reduced" : + return escript.Data() + elif name == "y_contact_reduced" : + return escript.Data() elif name == "r" : return escript.Data() elif name == "q" : @@ -1709,8 +1974,10 @@ self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR), "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR), "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), + "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR), "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), + "g_reduced": PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE), "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH), "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)} self.setSymmetryOn() @@ -1772,6 +2039,26 @@ return escript.Data() elif name == "y_contact" : return escript.Data() + elif name == "A_reduced" : + return escript.Data() + elif name == "B_reduced" : + return escript.Data() + elif name == "C_reduced" : + return escript.Data() + elif name == "D_reduced" : + return escript.Data() + elif name == "X_reduced" : + return escript.Data() + elif name == "Y_reduced" : + return self.getCoefficient("f_reduced") + elif name == "d_reduced" : + return escript.Data() + elif name == "y_reduced" : + return self.getCoefficient("g_reduced") + elif name == "d_contact_reduced" : + return escript.Data() + elif name == "y_contact_reduced" : + return escript.Data() elif name == "r" : return self.getCoefficient("r") elif name == "q" : @@ -1783,11 +2070,11 @@ """ Class to define a Lame equation problem: - M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[j])[i])[j] = F_i -grad(S{sigma}[i,j])[j] } + M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] } with natural boundary conditons: - M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) - S{lambda}*grad(u[j])[i]) = f_i -n[j]*S{sigma}[i,j] } + M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] } and constraints: @@ -1807,7 +2094,7 @@ "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)} self.setSymmetryOn() - def setValue(self,**coefficients): + def setValues(self,**coefficients): """ sets new values to coefficients @@ -1830,7 +2117,7 @@ depending of reduced order is used for the representation of the equation. @raise IllegalCoefficient: if an unknown coefficient keyword is used. """ - super(LameEquation, self).setValue(**coefficients) + super(LameEquation, self).setValues(**coefficients) def getCoefficientOfGeneralPDE(self,name): """ @@ -1870,492 +2157,25 @@ 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 - -class AdvectivePDE(LinearPDE): - """ - In cases of PDEs dominated by the advection terms M{B} and M{C} against the adevctive terms M{A} - up-winding has been used. The L{AdvectivePDE} class applies SUPG upwinding to the advective terms. - - In the following we set - - M{Z[j]=C[j]-B[j]} - - or - - M{Z[i,k,l]=C[i,k,l]-B[i,l,k]} - - To measure the dominance of the advective terms over the diffusive term M{A} the - X{Pelclet number} M{P} is used. It is defined as - - M{P=h|Z|/(2|A|)} - - where M{|.|} denotes the L{length} of the arument and M{h} is the local cell size - from L{getSize}. Where M{|A|==0} M{P} is M{S{infinity}}. - - From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated: - - M{S{Xi}=S{xi}(P) h/|Z|} - - where M{S{xi}} is a suitable function of the Peclet number. - - In the case of a single PDE the coefficient are up-dated in the following way: - - M{A[i,j] S{<-} A[i,j] + S{Xi} * Z[j] * Z[l]} - - M{B[j] S{<-} B[j] + S{Xi} * C[j] * D} - - M{C[j] S{<-} C[j] + S{Xi} * B[j] * D} - - M{X[j] S{<-} X[j] + S{Xi} * Z[j] * Y} - - Similar for the case of a systems of PDEs: - - M{A[i,j,k,l] S{<-} A[i,j,k,l]+ S{delta}[p,m] * S{Xi} * Z[p,i,j] * Z[m,k,l]} - - M{B[i,j,k] S{<-} B[i,j,k] + S{delta}[p,m] * S{Xi} * D[p,k] * C[m,i,j]} - - M{C[i,k,l] S{<-} C[i,k,l] + S{delta}[p,m] * S{Xi} * D[p,k] * B[m,l,i]} - - M{X[i,j] S{<-} X[i,j] + S{delta}[p,m] * S{Xi} * Y[p] * Z[m,i,j]} - - where M{S{delta}} is L{kronecker}. - Using upwinding in this form, introduces an additonal error which is proprtional to the cell size M{h} - but with the intension to stabilize the solution. - - """ - def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False): - """ - creates a linear, steady, second order PDE on a L{Domain} - - @param domain: domain of the PDE - @type domain: L{Domain} - @param numEquations: number of equations. If numEquations==None the number of equations - is exracted from the PDE coefficients. - @param numSolutions: number of solution components. If numSolutions==None the number of solution components - is exracted from the PDE coefficients. - @param xi: defines a function which returns for any given Preclet number as L{Scalar} object the - M{S{xi}}-value used to define the stabilization parameters. If equal to None, L{ELMAN_RAMAGE} is used. - @type xi: callable object which returns a L{Scalar} object. - @param debug: if True debug informations are printed. - """ - super(AdvectivePDE, self).__init__(domain,\ - numEquations,numSolutions,debug) - if xi==None: - self.__xi=AdvectivePDE.ELMAN_RAMAGE - else: - self.__xi=xi - self.__Xi=escript.Data() - - def setValue(**coefficients): - """ - sets new values to coefficients - - @param coefficients: new values assigned to coefficients - @keyword A: value for coefficient A. - @type A: any type that can be casted to L{Data} object on L{Function}. - @keyword B: value for coefficient B - @type B: any type that can be casted to L{Data} object on L{Function}. - @keyword C: value for coefficient C - @type C: any type that can be casted to L{Data} object on L{Function}. - @keyword D: value for coefficient D - @type D: any type that can be casted to L{Data} object on L{Function}. - @keyword X: value for coefficient X - @type X: any type that can be casted to L{Data} object on L{Function}. - @keyword Y: value for coefficient Y - @type Y: any type that can be casted to L{Data} object on L{Function}. - @keyword d: value for coefficient d - @type d: any type that can be casted to L{Data} object on L{FunctionOnBoundary}. - @keyword y: value for coefficient y - @type y: any type that can be casted to L{Data} object on L{FunctionOnBoundary}. - @keyword d_contact: value for coefficient d_contact - @type d_contact: any type that can be casted to L{Data} object on L{FunctionOnContactOne}. - or L{FunctionOnContactZero}. - @keyword y_contact: value for coefficient y_contact - @type y_contact: any type that can be casted to L{Data} object on L{FunctionOnContactOne}. - or L{FunctionOnContactZero}. - @keyword r: values prescribed to the solution at the locations of constraints - @type r: any type that can be casted to L{Data} object on L{Solution} or L{ReducedSolution} - depending of reduced order is used for the solution. - @keyword q: mask for location of constraints - @type q: any type that can be casted to L{Data} object on L{Solution} or L{ReducedSolution} - depending of reduced order is used for the representation of the equation. - @raise IllegalCoefficient: if an unknown coefficient keyword is used. - - """ - if "A" in coefficients.keys() or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data() - super(AdvectivePDE, self).setValue(**coefficients) - - def ELMAN_RAMAGE(self,P): - """ - Predefined function to set a values for M{S{xi}} from a Preclet number M{P}. - This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002) - - M{S{xi}(P)=0} for M{P<1} - - M{S{xi}(P)=(1-1/P)/2} otherwise - - @param P: Preclet number - @type P: L{Scalar} - @return: up-wind weightimg factor - @rtype: L{Scalar} - """ - return util.wherePositive(P-1.)*0.5*(1.-1./(P+1.e-15)) - - def SIMPLIFIED_BROOK_HUGHES(self,P): - """ - Predefined function to set a values for M{S{xi}} from a Preclet number M{P}. - The original methods is - - M{S{xi}(P)=coth(P)-1/P} - - As the evaluation of M{coth} is expensive we are using the approximation: - - - M{S{xi}(P)=P/3} where M{P<3} - - M{S{xi}(P)=1/2} otherwise - - @param P: Preclet number - @type P: L{Scalar} - @return: up-wind weightimg factor - @rtype: L{Scalar} - """ - c=util.whereNegative(P-3.) - return P/6.*c+1./2.*(1.-c) - - def HALF(self,P): - """ - Predefined function to set value M{1/2} for M{S{xi}} - - @param P: Preclet number - @type P: L{Scalar} - @return: up-wind weightimg factor - @rtype: L{Scalar} - """ - return escript.Scalar(0.5,P.getFunctionSpace()) - - 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. - - def __getXi(self): - if self.__Xi.isEmpty(): - B=self.getCoefficient("B") - C=self.getCoefficient("C") - A=self.getCoefficient("A") - h=self.getDomain().getSize() - self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A")) - if not C.isEmpty() or not B.isEmpty(): - if not C.isEmpty() and not B.isEmpty(): - flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A")) - if self.getNumEquations()>1: - if self.getNumSolutions()>1: - for i in range(self.getNumEquations()): - for k in range(self.getNumSolutions()): - for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2 - # flux=C-util.reorderComponents(B,[0,2,1]) - else: - for i in range(self.getNumEquations()): - for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2 - # flux=C-B - else: - if self.getNumSolutions()>1: - for k in range(self.getNumSolutions()): - for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2 - # flux=C-util.reorderComponents(B,[1,0]) - else: - for l in range(self.getDim()): flux2+=(C[l]-B[l])**2 - #flux=C-B - length_of_flux=util.sqrt(flux2) - elif C.isEmpty(): - length_of_flux=util.length(B) - #flux=B - else: - length_of_flux=util.length(C) - #flux=C - - #length_of_flux=util.length(flux) - flux_max=util.Lsup(length_of_flux) - if flux_max>0.: - # length_of_A=util.inner(flux,util.tensormutiply(A,flux)) - length_of_A=util.length(A) - A_max=util.Lsup(length_of_A) - if A_max>0: - inv_A=1./(length_of_A+A_max*self.__TOL) - else: - inv_A=1./self.__TOL - peclet_number=length_of_flux*h/2*inv_A - xi=self.__xi(peclet_number) - self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL) - self.trace("preclet number = %e"%util.Lsup(peclet_number)) - return self.__Xi - - - 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} - @raise IllegalCoefficient: if name is not one of coefficients - M{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 not self.getNumEquations() == self.getNumSolutions(): - raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal." - - if name == "A" : - A=self.getCoefficient("A") - B=self.getCoefficient("B") - C=self.getCoefficient("C") - if B.isEmpty() and C.isEmpty(): - Aout=A - else: - if A.isEmpty(): - Aout=self.createNewCoefficient("A") - else: - Aout=A[:] - Xi=self.__getXi() - if self.getNumEquations()>1: - for i in range(self.getNumEquations()): - for j in range(self.getDim()): - for k in range(self.getNumSolutions()): - for l in range(self.getDim()): - if not C.isEmpty() and not B.isEmpty(): - # tmp=C-util.reorderComponents(B,[0,2,1]) - # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1) - for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k]) - elif C.isEmpty(): - for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k] - # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1) - else: - for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l] - # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1) - else: - for j in range(self.getDim()): - for l in range(self.getDim()): - if not C.isEmpty() and not B.isEmpty(): - Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l]) - elif C.isEmpty(): - Aout[j,l]+=Xi*B[j]*B[l] - else: - 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) - return Aout - elif name == "B" : - B=self.getCoefficient("B") - C=self.getCoefficient("C") - D=self.getCoefficient("D") - if C.isEmpty() or D.isEmpty(): - Bout=B - else: - Xi=self.__getXi() - if B.isEmpty(): - Bout=self.createNewCoefficient("B") - else: - Bout=B[:] - if self.getNumEquations()>1: - for k in range(self.getNumSolutions()): - for p in range(self.getNumEquations()): - tmp=Xi*D[p,k] - for i in range(self.getNumEquations()): - for j in range(self.getDim()): - Bout[i,j,k]+=tmp*C[p,i,j] - # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1) - else: - tmp=Xi*D - for j in range(self.getDim()): Bout[j]+=tmp*C[j] - # Bout=Bout+Xi*D*C - return Bout - elif name == "C" : - B=self.getCoefficient("B") - C=self.getCoefficient("C") - D=self.getCoefficient("D") - if B.isEmpty() or D.isEmpty(): - Cout=C - else: - Xi=self.__getXi() - if C.isEmpty(): - Cout=self.createNewCoefficient("C") - else: - Cout=C[:] - if self.getNumEquations()>1: - for k in range(self.getNumSolutions()): - for p in range(self.getNumEquations()): - tmp=Xi*D[p,k] - for i in range(self.getNumEquations()): - for l in range(self.getDim()): - Cout[i,k,l]+=tmp*B[p,l,i] - # Cout=Cout+Xi*B[p,l,i]*D[p,k] - else: - tmp=Xi*D - for j in range(self.getDim()): Cout[j]+=tmp*B[j] - # Cout=Cout+tmp*D*B - return Cout - elif name == "D" : - return self.getCoefficient("D") - elif name == "X" : - X=self.getCoefficient("X") - Y=self.getCoefficient("Y") - B=self.getCoefficient("B") - C=self.getCoefficient("C") - if Y.isEmpty() or (B.isEmpty() and C.isEmpty()): - Xout=X - else: - if X.isEmpty(): - Xout=self.createNewCoefficient("X") - else: - Xout=X[:] - Xi=self.__getXi() - if self.getNumEquations()>1: - for p in range(self.getNumEquations()): - tmp=Xi*Y[p] - for i in range(self.getNumEquations()): - for j in range(self.getDim()): - if not C.isEmpty() and not B.isEmpty(): - Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i]) - # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1) - elif C.isEmpty(): - Xout[i,j]-=tmp*B[p,j,i] - # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1) - else: - Xout[i,j]+=tmp*C[p,i,j] - # Xout=X_out+Xi*util.inner(Y,C,offset=1) - else: - tmp=Xi*Y - for j in range(self.getDim()): - if not C.isEmpty() and not B.isEmpty(): - Xout[j]+=tmp*(C[j]-B[j]) - # Xout=Xout+Xi*Y*(C-B) - elif C.isEmpty(): - Xout[j]-=tmp*B[j] - # Xout=Xout-Xi*Y*B - else: - Xout[j]+=tmp*C[j] - # Xout=Xout+Xi*Y*C - return Xout - elif name == "Y" : - return self.getCoefficient("Y") - elif name == "d" : - return self.getCoefficient("d") - elif name == "y" : - return self.getCoefficient("y") - elif name == "d_contact" : - return self.getCoefficient("d_contact") - elif name == "y_contact" : - return self.getCoefficient("y_contact") - 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 - -class AdvectionDiffusion(LinearPDE): - """ - Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form - - M{S{omega}*u + inner(v,grad(u))- grad(matrixmult(k_bar,grad(u))[j])[j] = f} - - 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} - @param debug: if True debug informations are printed. - - """ - super(AdvectionDiffusion, self).__init__(domain,1,1,debug) - 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} object on L{Function}. - @keyword k: value for coefficient M{k} - @type k: any type that can be casted to L{Tensor} object on L{Function}. - @keyword v: value for coefficient M{v} - @type v: any type that can be casted to L{Vector} object on L{Function}. - @keyword upwind: value for upwind term M{upwind} - @type upwind: any type that can be casted to L{Vector} object on L{Function}. - @keyword f: value for right hand side M{f} - @type f: any type that can be casted to L{Scalar} object on L{Function}. - @keyword alpha: value for right hand side M{S{alpha}} - @type alpha: any type that can be casted to L{Scalar} object on L{FunctionOnBoundary}. - @keyword g: value for right hand side M{g} - @type g: any type that can be casted to L{Scalar} object on L{FunctionOnBoundary}. - @keyword r: prescribed values M{r} for the solution in constraints. - @type r: any type that can be casted to L{Scalar} object on L{Solution} or L{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} object on L{Solution} or L{ReducedSolution} - depending of reduced order is used for the representation of the equation. - @raise IllegalCoefficient: if an unknown coefficient keyword is used. - """ - super(AdvectionDiffusion, self).setValue(**coefficients) - - 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} - @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" : + elif name == "A_reduced" : return escript.Data() - elif name == "C" : - return self.getCoefficient("v") - elif name == "D" : - return self.getCoefficient("omega") - elif name == "X" : + elif name == "B_reduced" : 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" : + elif name == "C_reduced" : return escript.Data() - elif name == "y_contact" : + elif name == "D_reduced" : + return escript.Data() + elif name == "X_reduced" : + return escript.Data() + elif name == "Y_reduced" : + return escript.Data() + elif name == "d_reduced" : + return escript.Data() + elif name == "y_reduced" : + return escript.Data() + elif name == "d_contact_reduced" : + return escript.Data() + elif name == "y_contact_reduced" : return escript.Data() elif name == "r" : return self.getCoefficient("r") @@ -2364,189 +2184,3 @@ else: raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name - -# $Log$ -# Revision 1.14 2005/09/22 01:54:57 jgs -# Merge of development branch dev-02 back to main trunk on 2005-09-22 -# -# Revision 1.13 2005/09/15 03:44:19 jgs -# Merge of development branch dev-02 back to main trunk on 2005-09-15 -# -# Revision 1.12 2005/09/01 03:31:28 jgs -# Merge of development branch dev-02 back to main trunk on 2005-09-01 -# -# Revision 1.11 2005/08/23 01:24:28 jgs -# Merge of development branch dev-02 back to main trunk on 2005-08-23 -# -# Revision 1.10 2005/08/12 01:45:36 jgs -# erge of development branch dev-02 back to main trunk on 2005-08-12 -# -# Revision 1.9.2.17 2005/09/21 07:03:33 matt -# PDECoefficient and LinearPDE are now new style classes (introduced in Python -# 2.2). Classes Poisson, Helmholtz, LameEquation and AdvectivePDE have been -# modified to instead use portable/cooperative "super" calls to extend base -# class methods. -# -# Revision 1.9.2.16 2005/09/16 01:54:37 matt -# Removed redundant if-loop. -# -# Revision 1.9.2.15 2005/09/14 08:09:18 matt -# Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs. -# -# Revision 1.9.2.14 2005/09/07 06:26:16 gross -# the solver from finley are put into the standalone package paso now -# -# Revision 1.9.2.13 2005/08/31 08:45:03 gross -# in the case of lumping no new system is allocated if the constraint is changed. -# -# Revision 1.9.2.12 2005/08/31 07:10:23 gross -# test for Lumping added -# -# Revision 1.9.2.11 2005/08/30 01:53:45 gross -# bug in format fixed. -# -# Revision 1.9.2.10 2005/08/26 07:14:17 gross -# a few more bugs in linearPDE fixed. remaining problem are finley problems -# -# Revision 1.9.2.9 2005/08/26 06:30:45 gross -# fix for reported bug 0000004. test_linearPDE passes a few more tests -# -# Revision 1.9.2.8 2005/08/26 04:30:13 gross -# gneric unit testing for linearPDE -# -# Revision 1.9.2.7 2005/08/25 07:06:50 gross -# linearPDE documentation is parsed now by epydoc. there is still a problem with links into escriptcpp.so -# -# Revision 1.9.2.6 2005/08/24 05:01:24 gross -# problem with resetting the matrix in case of resetting its values to 0 fixed. -# -# Revision 1.9.2.5 2005/08/24 02:03:28 gross -# epydoc mark up partially fixed -# -# Revision 1.9.2.4 2005/08/22 07:11:09 gross -# some problems with LinearPDEs fixed. -# -# Revision 1.9.2.3 2005/08/18 04:48:48 gross -# the methods SetLumping*() are removed. Lumping is set trough setSolverMethod(LinearPDE.LUMPING) -# -# Revision 1.9.2.2 2005/08/18 04:39:32 gross -# the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now -# -# Revision 1.9.2.1 2005/07/29 07:10:27 gross -# new functions in util and a new pde type in linearPDEs -# -# Revision 1.1.2.25 2005/07/28 04:21:09 gross -# Lame equation: (linear elastic, isotropic) added -# -# Revision 1.1.2.24 2005/07/22 06:37:11 gross -# some extensions to modellib and linearPDEs -# -# Revision 1.1.2.23 2005/05/13 00:55:20 cochrane -# Fixed up some docstrings. Moved module-level functions to top of file so -# that epydoc and doxygen can pick them up properly. -# -# Revision 1.1.2.22 2005/05/12 11:41:30 gross -# some basic Models have been added -# -# Revision 1.1.2.21 2005/05/12 07:16:12 cochrane -# Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of -# file so that the AdvectivePDE class is picked up by doxygen. Some -# reformatting of docstrings. Addition of code to make equations come out -# as proper LaTeX. -# -# Revision 1.1.2.20 2005/04/15 07:09:08 gross -# some problems with functionspace and linearPDEs fixed. -# -# Revision 1.1.2.19 2005/03/04 05:27:07 gross -# bug in SystemPattern fixed. -# -# Revision 1.1.2.18 2005/02/08 06:16:45 gross -# Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed -# -# Revision 1.1.2.17 2005/02/08 05:56:19 gross -# Reference Number handling added -# -# Revision 1.1.2.16 2005/02/07 04:41:28 gross -# some function exposed to python to make mesh merging running -# -# Revision 1.1.2.15 2005/02/03 00:14:44 gross -# timeseries add and ESySParameter.py renames esysXML.py for consistence -# -# Revision 1.1.2.14 2005/02/01 06:44:10 gross -# new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working -# -# Revision 1.1.2.13 2005/01/25 00:47:07 gross -# updates in the documentation -# -# Revision 1.1.2.12 2005/01/12 01:28:04 matt -# Added createCoefficient method for linearPDEs. -# -# Revision 1.1.2.11 2005/01/11 01:55:34 gross -# a problem in linearPDE class fixed -# -# Revision 1.1.2.10 2005/01/07 01:13:29 gross -# some bugs in linearPDE fixed -# -# Revision 1.1.2.9 2005/01/06 06:24:58 gross -# some bugs in slicing fixed -# -# Revision 1.1.2.8 2005/01/05 04:21:40 gross -# FunctionSpace checking/matchig in slicing added -# -# Revision 1.1.2.7 2004/12/29 10:03:41 gross -# bug in setValue fixed -# -# Revision 1.1.2.6 2004/12/29 05:29:59 gross -# AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data() -# -# Revision 1.1.2.5 2004/12/29 00:18:41 gross -# AdvectivePDE added -# -# Revision 1.1.2.4 2004/12/24 06:05:41 gross -# some changes in linearPDEs to add AdevectivePDE -# -# Revision 1.1.2.3 2004/12/16 00:12:34 gross -# __init__ of LinearPDE does not accept any coefficient anymore -# -# Revision 1.1.2.2 2004/12/14 03:55:01 jgs -# *** empty log message *** -# -# Revision 1.1.2.1 2004/12/12 22:53:47 gross -# linearPDE has been renamed LinearPDE -# -# Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross -# GMRES added -# -# Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross -# options for GMRES and PRES20 added -# -# Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross -# some small changes -# -# Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross -# Finley solves 4M unknowns now -# -# Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross -# poisson solver added -# -# Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross -# a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry -# -# Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross -# finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed -# -# Revision 1.1.1.1 2004/10/26 06:53:56 jgs -# initial import of project esys2 -# -# Revision 1.3.2.3 2004/10/26 06:43:48 jgs -# committing Lutz's and Paul's changes to brach jgs -# -# Revision 1.3.4.1 2004/10/20 05:32:51 cochrane -# Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form. -# -# Revision 1.3 2004/09/23 00:53:23 jgs -# minor fixes -# -# Revision 1.1 2004/08/28 12:58:06 gross -# SimpleSolve is not running yet: problem with == of functionsspace -#