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

revision 3521 by lgao, Fri May 13 04:15:29 2011 UTC revision 3522 by gross, Tue May 24 00:57:58 2011 UTC
# Line 1243  class PDECoef(object): Line 1243  class PDECoef(object):
1243                      the PDE                      the PDE
1244      :cvar REDUCED: indicator that coefficient is defined through a reduced      :cvar REDUCED: indicator that coefficient is defined through a reduced
1245                     solution of the PDE                     solution of the PDE
1246        :cvar DIRACDELTA: indicator that coefficient is defined as Dirac delta functions
1247      :cvar BY_EQUATION: indicator that the dimension of the coefficient shape is      :cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
1248                         defined by the number of PDE equations                         defined by the number of PDE equations
1249      :cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is      :cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
# Line 1271  class PDECoef(object): Line 1272  class PDECoef(object):
1272      INTERIOR_REDUCED=13      INTERIOR_REDUCED=13
1273      BOUNDARY_REDUCED=14      BOUNDARY_REDUCED=14
1274      CONTACT_REDUCED=15      CONTACT_REDUCED=15
1275        DIRACDELTA=16
1276
1277      def __init__(self, where, pattern, altering):      def __init__(self, where, pattern, altering):
1278         """         """
# Line 1279  class PDECoef(object): Line 1281  class PDECoef(object):
1281         :param where: describes where the coefficient lives         :param where: describes where the coefficient lives
1282         :type where: one of `INTERIOR`, `BOUNDARY`, `CONTACT`, `SOLUTION`,         :type where: one of `INTERIOR`, `BOUNDARY`, `CONTACT`, `SOLUTION`,
1283                      `REDUCED`, `INTERIOR_REDUCED`, `BOUNDARY_REDUCED`,                      `REDUCED`, `INTERIOR_REDUCED`, `BOUNDARY_REDUCED`,
1284                      `CONTACT_REDUCED`                      `CONTACT_REDUCED`, 'DIRACDELTA'
1285         :param pattern: describes the shape of the coefficient and how the shape         :param pattern: describes the shape of the coefficient and how the shape
1286                         is built for a given spatial dimension and numbers of                         is built for a given spatial dimension and numbers of
1287                         equations and solutions in then PDE. For instance,                         equations and solutions in then PDE. For instance,
# Line 1334  class PDECoef(object): Line 1336  class PDECoef(object):
1336              return escript.FunctionOnContactZero(domain)              return escript.FunctionOnContactZero(domain)
1337         elif self.what==self.CONTACT_REDUCED:         elif self.what==self.CONTACT_REDUCED:
1338              return escript.ReducedFunctionOnContactZero(domain)              return escript.ReducedFunctionOnContactZero(domain)
1339           elif self.what==self.DIRACDELTA:
1340                return escript.DiracDeltaFunctions(domain)
1341         elif self.what==self.SOLUTION:         elif self.what==self.SOLUTION:
1342              if reducedEquationOrder and reducedSolutionOrder:              if reducedEquationOrder and reducedSolutionOrder:
1343                  return escript.ReducedSolution(domain)                  return escript.ReducedSolution(domain)
# Line 2661  class LinearPDE(LinearProblem): Line 2665  class LinearPDE(LinearProblem):
2665         y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),         y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2666         d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),         d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2667         y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),         y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2668           d_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2669           y_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2670         r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),         r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2671         q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )         q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2672
# Line 2706  class LinearPDE(LinearProblem): Line 2712  class LinearPDE(LinearProblem):
2712        out=out and self.checkSymmetricTensor("d_reduced", verbose)        out=out and self.checkSymmetricTensor("d_reduced", verbose)
2713        out=out and self.checkSymmetricTensor("d_contact", verbose)        out=out and self.checkSymmetricTensor("d_contact", verbose)
2714        out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)        out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2715          out=out and self.checkSymmetricTensor("d_dirac", verbose)
2716        return out        return out
2717
2718     def createOperator(self):     def createOperator(self):
# Line 2774  class LinearPDE(LinearProblem): Line 2781  class LinearPDE(LinearProblem):
2781                   d=self.getCoefficient("d")                   d=self.getCoefficient("d")
2782                   D_reduced=self.getCoefficient("D_reduced")                   D_reduced=self.getCoefficient("D_reduced")
2783                   d_reduced=self.getCoefficient("d_reduced")                   d_reduced=self.getCoefficient("d_reduced")
2784                     d_dirac=self.getCoefficient("d_dirac")
2785
2786                   if not D.isEmpty():                   if not D.isEmpty():
2787                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
2788                          D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))                          D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))
# Line 2796  class LinearPDE(LinearProblem): Line 2805  class LinearPDE(LinearProblem):
2805                          D_reduced_times_e=D_reduced                          D_reduced_times_e=D_reduced
2806                   else:                   else:
2807                      D_reduced_times_e=escript.Data()                      D_reduced_times_e=escript.Data()
2808
2809                   if not d_reduced.isEmpty():                   if not d_reduced.isEmpty():
2810                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
2811                          d_reduced_times_e=util.matrix_mult(d_reduced,numpy.ones((self.getNumSolutions(),)))                          d_reduced_times_e=util.matrix_mult(d_reduced,numpy.ones((self.getNumSolutions(),)))
# Line 2803  class LinearPDE(LinearProblem): Line 2813  class LinearPDE(LinearProblem):
2813                          d_reduced_times_e=d_reduced                          d_reduced_times_e=d_reduced
2814                   else:                   else:
2815                      d_reduced_times_e=escript.Data()                      d_reduced_times_e=escript.Data()
2816
2817                     if not d_dirac.isEmpty():
2818                         if self.getNumSolutions()>1:
2819                            d_dirac_times_e=util.matrix_mult(d_dirac,numpy.ones((self.getNumSolutions(),)))
2820                         else:
2821                            d_reduced_dirac_e=d_dirac
2822                     else:
2823                        d_dirac_times_e=escript.Data()
2824
2825                   self.resetOperator()                   self.resetOperator()
2826                   operator=self.getCurrentOperator()                   operator=self.getCurrentOperator()
2828              hrz_lumping=( self.getSolverOptions().getSolverMethod() ==  SolverOptions.HRZ_LUMPING )              hrz_lumping=( self.getSolverOptions().getSolverMethod() ==  SolverOptions.HRZ_LUMPING )
2829              self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e,  hrz_lumping )              self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e, d_dirac_times_e,  hrz_lumping )
2831                   else:                   else:
2833                                                   escript.Data(), \                                                   escript.Data(), \
2834                                                   D_times_e, \                                                   D_times_e, \
2835                                                   d_times_e,\                                                   d_times_e,\
2836                                                   escript.Data())                                                   escript.Data(),\
2837                                                     d_dirac_times_e)
2839                                                   escript.Data(), \                                                   escript.Data(), \
2840                                                   D_reduced_times_e, \                                                   D_reduced_times_e, \
2841                                                   d_reduced_times_e,\                                                   d_reduced_times_e,\
2842                                                     escript.Data(), \
2843                                                   escript.Data())                                                   escript.Data())
2844                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
2845                if not self.isRightHandSideValid():                if not self.isRightHandSideValid():
# Line 2829  class LinearPDE(LinearProblem): Line 2849  class LinearPDE(LinearProblem):
2849                                 self.getCoefficient("X"), \                                 self.getCoefficient("X"), \
2850                                 self.getCoefficient("Y"),\                                 self.getCoefficient("Y"),\
2851                                 self.getCoefficient("y"),\                                 self.getCoefficient("y"),\
2852                                 self.getCoefficient("y_contact"))                                 self.getCoefficient("y_contact"), \
2853                                   self.getCoefficient("y_dirac"))
2855                                 self.getCoefficient("X_reduced"), \                                 self.getCoefficient("X_reduced"), \
2856                                 self.getCoefficient("Y_reduced"),\                                 self.getCoefficient("Y_reduced"),\
2857                                 self.getCoefficient("y_reduced"),\                                 self.getCoefficient("y_reduced"),\
2858                                 self.getCoefficient("y_contact_reduced"))                                 self.getCoefficient("y_contact_reduced"), \
2859                                   escript.Data())
2860                   self.trace("New right hand side has been built.")                   self.trace("New right hand side has been built.")
2861                   self.validRightHandSide()                   self.validRightHandSide()
2862                self.insertConstraint(rhs_only=False)                self.insertConstraint(rhs_only=False)
# Line 2855  class LinearPDE(LinearProblem): Line 2877  class LinearPDE(LinearProblem):
2877                                 self.getCoefficient("d"), \                                 self.getCoefficient("d"), \
2878                                 self.getCoefficient("y"), \                                 self.getCoefficient("y"), \
2879                                 self.getCoefficient("d_contact"), \                                 self.getCoefficient("d_contact"), \
2880                                 self.getCoefficient("y_contact"))                                 self.getCoefficient("y_contact"), \
2881                                   self.getCoefficient("d_dirac"), \
2882                                   self.getCoefficient("y_dirac"))
2884                                 self.getCoefficient("A_reduced"), \                                 self.getCoefficient("A_reduced"), \
2885                                 self.getCoefficient("B_reduced"), \                                 self.getCoefficient("B_reduced"), \
# Line 2866  class LinearPDE(LinearProblem): Line 2890  class LinearPDE(LinearProblem):
2890                                 self.getCoefficient("d_reduced"), \                                 self.getCoefficient("d_reduced"), \
2891                                 self.getCoefficient("y_reduced"), \                                 self.getCoefficient("y_reduced"), \
2892                                 self.getCoefficient("d_contact_reduced"), \                                 self.getCoefficient("d_contact_reduced"), \
2893                                 self.getCoefficient("y_contact_reduced"))                                 self.getCoefficient("y_contact_reduced"), \
2894                                   escript.Data(), \
2895                                   escript.Data())
2896                   self.insertConstraint(rhs_only=False)                   self.insertConstraint(rhs_only=False)
2897                   self.trace("New system has been built.")                   self.trace("New system has been built.")
2898                   self.validOperator()                   self.validOperator()
# Line 2878  class LinearPDE(LinearProblem): Line 2904  class LinearPDE(LinearProblem):
2904                                 self.getCoefficient("X"), \                                 self.getCoefficient("X"), \
2905                                 self.getCoefficient("Y"),\                                 self.getCoefficient("Y"),\
2906                                 self.getCoefficient("y"),\                                 self.getCoefficient("y"),\
2907                                 self.getCoefficient("y_contact"))                                 self.getCoefficient("y_contact"), \
2908                                   self.getCoefficient("y_dirac") )
2910                                 self.getCoefficient("X_reduced"), \                                 self.getCoefficient("X_reduced"), \
2911                                 self.getCoefficient("Y_reduced"),\                                 self.getCoefficient("Y_reduced"),\
2912                                 self.getCoefficient("y_reduced"),\                                 self.getCoefficient("y_reduced"),\
2913                                 self.getCoefficient("y_contact_reduced"))                                 self.getCoefficient("y_contact_reduced"), \
2914                                   escript.Data())
2915                   self.insertConstraint(rhs_only=True)                   self.insertConstraint(rhs_only=True)
2916                   self.trace("New right hand side has been built.")                   self.trace("New right hand side has been built.")
2917                   self.validRightHandSide()                   self.validRightHandSide()
# Line 2900  class LinearPDE(LinearProblem): Line 2928  class LinearPDE(LinearProblem):
2928                              self.getCoefficient("d"), \                              self.getCoefficient("d"), \
2929                              escript.Data(),\                              escript.Data(),\
2930                              self.getCoefficient("d_contact"), \                              self.getCoefficient("d_contact"), \
2931                                escript.Data(),                   \
2932                                self.getCoefficient("d_dirac"),   \
2933                              escript.Data())                              escript.Data())
2935                              self.getCoefficient("A_reduced"), \                              self.getCoefficient("A_reduced"), \
# Line 2911  class LinearPDE(LinearProblem): Line 2941  class LinearPDE(LinearProblem):
2941                              self.getCoefficient("d_reduced"), \                              self.getCoefficient("d_reduced"), \
2942                              escript.Data(),\                              escript.Data(),\
2943                              self.getCoefficient("d_contact_reduced"), \                              self.getCoefficient("d_contact_reduced"), \
2944                                escript.Data(),  \
2945                                escript.Data(),  \
2946                              escript.Data())                              escript.Data())
2947                   self.insertConstraint(rhs_only=False)                   self.insertConstraint(rhs_only=False)
2948                   self.trace("New operator has been built.")                   self.trace("New operator has been built.")
# Line 3015  class LinearPDE(LinearProblem): Line 3047  class LinearPDE(LinearProblem):
3047        :type y_contact_reduced: any type that can be cast to a `Data`        :type y_contact_reduced: any type that can be cast to a `Data`
3048                                 object on `ReducedFunctionOnContactOne`                                 object on `ReducedFunctionOnContactOne`
3049                                 or `ReducedFunctionOnContactZero`                                 or `ReducedFunctionOnContactZero`
3050          :keyword d_dirac: value for coefficient ``d_dirac``
3051          :type d_dirac: any type that can be cast to a `Data` object on `DiracDeltaFunctions`
3052          :keyword y_dirac: value for coefficient ``y_dirac``
3053          :type y_dirac: any type that can be cast to a `Data` object on `DiracDeltaFunctions`
3054        :keyword r: values prescribed to the solution at the locations of        :keyword r: values prescribed to the solution at the locations of
3055                    constraints                    constraints
3056        :type r: any type that can be cast to a `Data` object on        :type r: any type that can be cast to a `Data` object on
# Line 3473  class TransportPDE(LinearProblem): Line 3509  class TransportPDE(LinearProblem):
3509        - *m_reduced[i,k]=m_reduced[k,i]*        - *m_reduced[i,k]=m_reduced[k,i]*
3510        - *d[i,k]=d[k,i]*        - *d[i,k]=d[k,i]*
3511        - *d_reduced[i,k]=d_reduced[k,i]*        - *d_reduced[i,k]=d_reduced[k,i]*
3512          - *d_dirac[i,k]=d_dirac[k,i]*
3513
3514     `TransportPDE` also supports solution discontinuities over a contact region     `TransportPDE` also supports solution discontinuities over a contact region
3515     in the domain. To specify the conditions across the discontinuity we are     in the domain. To specify the conditions across the discontinuity we are
# Line 3573  class TransportPDE(LinearProblem): Line 3610  class TransportPDE(LinearProblem):
3610         y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),         y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3611         d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),         d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3612         y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),         y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3613           d_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
3614           y_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3615         r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),         r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
3616         q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )         q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
3617
# Line 3620  class TransportPDE(LinearProblem): Line 3659  class TransportPDE(LinearProblem):
3659        out=out and self.checkSymmetricTensor("d_reduced", verbose)        out=out and self.checkSymmetricTensor("d_reduced", verbose)
3660        out=out and self.checkSymmetricTensor("d_contact", verbose)        out=out and self.checkSymmetricTensor("d_contact", verbose)
3661        out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)        out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
3662          out=out and self.checkSymmetricTensor("d_dirac", verbose)
3663        return out        return out
3664
3665     def setValue(self,**coefficients):     def setValue(self,**coefficients):
# Line 3694  class TransportPDE(LinearProblem): Line 3734  class TransportPDE(LinearProblem):
3734                         object on `FunctionOnContactOne` or `FunctionOnContactZero`                         object on `FunctionOnContactOne` or `FunctionOnContactZero`
3735        :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``        :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``
3736        :type y_contact_reduced: any type that can be cast to a `Data` object on `ReducedFunctionOnContactOne` or `ReducedFunctionOnContactZero`        :type y_contact_reduced: any type that can be cast to a `Data` object on `ReducedFunctionOnContactOne` or `ReducedFunctionOnContactZero`
3737
3738          :keyword d_dirac: value for coefficient ``d_dirac``
3739          :type d_dirac: any type that can be cast to a `Data` object on `DiracDeltaFunctions`
3740          :keyword y_dirac: value for coefficient ``y_dirac``
3741          :type y_dirac: any type that can be cast to a `Data` object on `DiracDeltaFunctions`
3742
3743        :keyword r: values prescribed to the solution at the locations of constraints        :keyword r: values prescribed to the solution at the locations of constraints
3744        :type r: any type that can be cast to a `Data` object on        :type r: any type that can be cast to a `Data` object on
3745                 `Solution` or `ReducedSolution`                 `Solution` or `ReducedSolution`
# Line 3849  class TransportPDE(LinearProblem): Line 3895  class TransportPDE(LinearProblem):
3895                              self.getCoefficient("d"),                              self.getCoefficient("d"),
3896                              self.getCoefficient("y"),                              self.getCoefficient("y"),
3897                              self.getCoefficient("d_contact"),                              self.getCoefficient("d_contact"),
3898                              self.getCoefficient("y_contact"))                              self.getCoefficient("y_contact"),
3899                                self.getCoefficient("d_dirac"),
3900                                self.getCoefficient("y_dirac") )
3902                              operator,                              operator,
3903                              righthandside,                              righthandside,
# Line 3863  class TransportPDE(LinearProblem): Line 3911  class TransportPDE(LinearProblem):
3911                              self.getCoefficient("d_reduced"),                              self.getCoefficient("d_reduced"),
3912                              self.getCoefficient("y_reduced"),                              self.getCoefficient("y_reduced"),
3913                              self.getCoefficient("d_contact_reduced"),                              self.getCoefficient("d_contact_reduced"),
3914                              self.getCoefficient("y_contact_reduced"))                              self.getCoefficient("y_contact_reduced"),
3915                                escript.Data(),
3916                                escript.Data() )
3917            operator.insertConstraint(righthandside,self.getCoefficient("q"),self.getCoefficient("r"),self.getConstraintWeightingFactor())            operator.insertConstraint(righthandside,self.getCoefficient("q"),self.getCoefficient("r"),self.getConstraintWeightingFactor())
3918            self.trace("New system has been built.")            self.trace("New system has been built.")
3919            self.validOperator()            self.validOperator()

Legend:
 Removed from v.3521 changed lines Added in v.3522