Sundance 2.1

00001 /* @HEADER@ */
00002 // ************************************************************************
00003 // 
00004 //                              Sundance
00005 //                 Copyright (2005) Sandia Corporation
00006 // 
00007 // Copyright (year first published) Sandia Corporation.  Under the terms 
00008 // of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 
00009 // retains certain rights in this software.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //                                                                                 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA                                                                                
00025 // Questions? Contact Kevin Long (krlong@sandia.gov), 
00026 // Sandia National Laboratories, Livermore, California, USA
00027 // 
00028 // ************************************************************************
00029 /* @HEADER@ */
00030 
00031 #include "SundanceExpr.hpp"
00032 #include "SundanceIntegral.hpp"
00033 #include "SundanceStdSumTransformations.hpp"
00034 
00035 #include "SundanceSumExpr.hpp"
00036 #include "SundanceProductExpr.hpp"
00037 #include "SundanceConstantExpr.hpp"
00038 #include "SundanceCoordExpr.hpp"
00039 #include "SundanceDerivative.hpp"
00040 #include "SundanceDiffOp.hpp"
00041 
00042 
00043 #include "SundanceUnaryMinus.hpp"
00044 #include "SundanceZeroExpr.hpp"
00045 #include "SundanceFunctionalPolynomial.hpp"
00046 #include "SundanceSumOfIntegrals.hpp"
00047 #include "SundanceSumOfBCs.hpp"
00048 #include "SundanceNullCellFilterStub.hpp"
00049 #include "SundanceOut.hpp"
00050 #include "Teuchos_TimeMonitor.hpp"
00051 
00052 using namespace Sundance;
00053 using namespace Sundance;
00054 
00055 using namespace Teuchos;
00056 using namespace Sundance;
00057 
00058 static Time& polysumTimer() 
00059 {
00060   static RCP<Time> rtn 
00061     = TimeMonitor::getNewTimer("IdentifyPolynomialSum"); 
00062   return *rtn;
00063 }
00064 static Time& reorderSumTimer() 
00065 {
00066   static RCP<Time> rtn 
00067     = TimeMonitor::getNewTimer("ReorderSum"); 
00068   return *rtn;
00069 }
00070 
00071 static Time& removeUnaryMinusTimer() 
00072 {
00073   static RCP<Time> rtn 
00074     = TimeMonitor::getNewTimer("RemoveUnaryMinusFromSum"); 
00075   return *rtn;
00076 }
00077 
00078 static Time& removeZeroTimer() 
00079 {
00080   static RCP<Time> rtn 
00081     = TimeMonitor::getNewTimer("RemoveZeroFromSum"); 
00082   return *rtn;
00083 }
00084 
00085 static Time& sumConstantsTimer() 
00086 {
00087   static RCP<Time> rtn 
00088     = TimeMonitor::getNewTimer("SumConstants"); 
00089   return *rtn;
00090 }
00091 
00092 static Time& moveConstantsTimer() 
00093 {
00094   static RCP<Time> rtn 
00095     = TimeMonitor::getNewTimer("MoveConstantsToLeftOfSum"); 
00096   return *rtn;
00097 }
00098 
00099 static Time& rearrangeRightSumWithConstantTimer() 
00100 {
00101   static RCP<Time> rtn 
00102     = TimeMonitor::getNewTimer("RearrangeRightSumWithConstant"); 
00103   return *rtn;
00104 }
00105 
00106 static Time& rearrangeLeftSumWithConstantTimer() 
00107 {
00108   static RCP<Time> rtn 
00109     = TimeMonitor::getNewTimer("RearrangeLeftSumWithConstant"); 
00110   return *rtn;
00111 }
00112 
00113 static Time& sumIntegralsTimer() 
00114 {
00115   static RCP<Time> rtn 
00116     = TimeMonitor::getNewTimer("SumIntegrals"); 
00117   return *rtn;
00118 }
00119 
00120 
00121 
00122 
00123 StdSumTransformations::StdSumTransformations()
00124   : SumTransformationSequence()
00125 {
00126   append(rcp(new RemoveZeroFromSum()));
00127 //  append(rcp(new RemoveUnaryMinusFromSum()));
00128   //  append(rcp(new ReorderSum()));
00129 
00130   append(rcp(new SumConstants()));
00131 
00132     append(rcp(new MoveConstantsToLeftOfSum()));
00133     append(rcp(new RearrangeRightSumWithConstant()));
00134     append(rcp(new RearrangeLeftSumWithConstant()));
00135   append(rcp(new IdentifyPolynomialSum()));
00136   append(rcp(new SumIntegrals()));
00137 }
00138 
00139 
00140 bool IdentifyPolynomialSum::doTransform(const RCP<ScalarExpr>& left, 
00141                                         const RCP<ScalarExpr>& right,
00142                                         int sign, 
00143                                         RCP<ScalarExpr>& rtn) const
00144 {
00145   TimeMonitor timer(polysumTimer());
00146   if (useOptimizedPolynomials())
00147     {
00148       /* see if this sum can be written as a polynomial in 
00149        * symbolic functions */
00150       if (FunctionalPolynomial::isConvertibleToPoly(left.get())
00151           && FunctionalPolynomial::isConvertibleToPoly(right.get()))
00152         {
00153           RCP<FunctionalPolynomial> lp = FunctionalPolynomial::toPoly(left);
00154           RCP<FunctionalPolynomial> rp = FunctionalPolynomial::toPoly(right);
00155           rtn = lp->addPoly(rp.get(), sign);
00156           return true;
00157         }
00158     }
00159   /* otherwise, do no transformation */
00160   return false;
00161 }
00162 
00163 
00164 bool ReorderSum::doTransform(const RCP<ScalarExpr>& left, 
00165                              const RCP<ScalarExpr>& right,
00166                              int sign, RCP<ScalarExpr>& rtn) const
00167 {
00168   TimeMonitor timer(reorderSumTimer());
00169   /* first we check to see whether the terms are already in order.
00170    * The left and right trees are already ordered, so that if the
00171    * first term on the right is after the last term on the left, the
00172    * combination of both is already ordered. In that case, nothing more needs
00173    * to be done. 
00174    */
00175   Expr L = Expr::handle(left);
00176   Expr R = Expr::handle(right);
00177   Sundance::Map<Expr, int> tree = L.getSumTree();
00178   Sundance::Map<Expr, int>::const_reverse_iterator iL = tree.rbegin();
00179   Expr endOfLeft = iL->first;
00180   Sundance::Map<Expr, int> rightTree = R.getSumTree();
00181   Sundance::Map<Expr, int>::const_iterator iR = rightTree.begin();
00182   Expr startOfRight = iR->first;
00183 
00184   if (endOfLeft.lessThan(startOfRight))
00185     {
00186       Tabs tab1;
00187       SUNDANCE_VERB_MEDIUM(tab1 << "Terms are already ordered, doing nothing");
00188       return false;
00189     }
00190   else
00191     {
00192       Tabs tab1;
00193 
00194       for (Map<Expr, int>::const_iterator i=rightTree.begin(); i!=rightTree.end(); i++)
00195         {
00196           int leftCount = 0;
00197           if (tree.containsKey(i->first))
00198             {
00199               leftCount = tree[i->first];
00200             }
00201           int count = leftCount + sign * i->second;
00202           tree.put(i->first, count);
00203         }
00204       SUNDANCE_VERB_MEDIUM(tab1 << "Ordered terms are: " << tree);      
00205 
00206 
00207       Expr sum = new ZeroExpr();
00208 
00209       /* add up all terms. If no terms are left after cancellations, the result will 
00210        * be zero. */
00211       for (Map<Expr, int>::const_iterator i=tree.begin(); i!=tree.end(); i++)
00212         {
00213           const Expr& term = i->first;
00214           int count = i->second;
00215           if (count==0) continue;
00216           if (count==1) sum = sum + term;
00217           else if (count==-1) sum = sum - term;
00218           else sum = sum + count*term;
00219         }
00220       
00221       rtn = getScalar(sum);
00222       return true;
00223     }
00224 }
00225 
00226 #ifdef OLD_CODE
00227 
00228 bool ReorderSum::doTransform(const RCP<ScalarExpr>& left, 
00229                              const RCP<ScalarExpr>& right,
00230                              int sign, RCP<ScalarExpr>& rtn) const
00231 {
00232   Tabs tabs;
00233   SUNDANCE_VERB_LOW(tabs << "trying ReorderSum");
00234 
00235   Tabs tab0;
00236   SUNDANCE_VERB_MEDIUM(tab0 << "L=" << l->toString());
00237   SUNDANCE_VERB_MEDIUM(tab0 << "R=" << r->toString());
00238   const SumExpr* sLeft = dynamic_cast<const SumExpr*>(l.get());
00239   const SumExpr* sRight = dynamic_cast<const SumExpr*>(r.get());
00240 
00241   if (sLeft==0 && sRight==0)
00242     {
00243       Tabs tab1;
00244       if (Expr::handle(r).lessThan(Expr::handle(l)))
00245         {
00246           if (sign>0)
00247             {
00248               SUNDANCE_VERB_MEDIUM(tab1 << "Both are non-sums: R < L, so reordering to R+L");
00249               rtn = getScalar(Expr::handle(r) + Expr::handle(l));
00250             }
00251           else
00252             {
00253               SUNDANCE_VERB_MEDIUM(tab1 << "Both are non-sums: R < L, reordering to -R+L");
00254               rtn = getScalar(-Expr::handle(r) + Expr::handle(l));
00255             }
00256           return true;
00257         }
00258     }
00259 
00260 
00261   if (sLeft != 0)
00262     {
00263       Tabs tab1;
00264       Expr ll = sLeft->left();
00265       Expr lr = sLeft->right();
00266       int lSign = sLeft->sign();
00267       if (Expr::handle(r).lessThan(ll))
00268         {
00269           if (sign > 0)
00270             {
00271               SUNDANCE_VERB_MEDIUM(tab1 << "L is a sum: R < LL, reordering to R+L");
00272               rtn = getScalar(Expr::handle(r) + Expr::handle(l));
00273             }
00274           else
00275             {
00276               SUNDANCE_VERB_MEDIUM(tab1 << "L is a sum: R < LL, reordering to -R+L");
00277               rtn = getScalar(-Expr::handle(r) + Expr::handle(l));
00278             }
00279           return true;
00280         }
00281       if (Expr::handle(r).lessThan(lr)) 
00282         {
00283           if (sign > 0)
00284             {
00285               if (lSign > 0)
00286                 {
00287                   SUNDANCE_VERB_MEDIUM(tab1 << "L is a sum: LL < R < LR, reordering to LL + R + LR");
00288                   rtn = getScalar(ll + Expr::handle(r) + lr);
00289                 }
00290               else
00291                 {
00292                   SUNDANCE_VERB_MEDIUM(tab1 << "L is a sum: LL < R < LR, reordering to LL + R - LR");
00293                   rtn = getScalar(ll + Expr::handle(r) - lr);
00294                 }
00295             }
00296           else
00297             {
00298               if (lSign > 0)
00299                 {
00300                   SUNDANCE_VERB_MEDIUM(tab1 << "L is a sum: LL < R < LR, reordering to LL - R + LR");
00301                   rtn = getScalar(ll - Expr::handle(r) + lr);
00302                 }
00303               else
00304                 {
00305                   SUNDANCE_VERB_MEDIUM(tab1 << "L is a sum: LL < R < LR, "
00306                                        "reordering to LL - R - LR");
00307                   rtn = getScalar(ll - Expr::handle(r) - lr);
00308                 }
00309             }
00310           return true;
00311         }
00312     }
00313   if (sRight != 0)
00314     {
00315       Tabs tab1;
00316       Expr rl = sRight->left();
00317       Expr rr = sRight->right();
00318       int rSign = sRight->sign();
00319 
00320       if (rr.lessThan(Expr::handle(l)))
00321         {
00322           if (sign > 0)
00323             {
00324               SUNDANCE_VERB_MEDIUM(tab1 << "R is a sum: R < L, reordering to R+L");
00325               rtn = getScalar(Expr::handle(r) + Expr::handle(l));
00326             }
00327           else
00328             {
00329               SUNDANCE_VERB_MEDIUM(tab1 << "R is a sum: R < L, reordering to -R+L");
00330               rtn = getScalar(-Expr::handle(r) + Expr::handle(l));
00331             }
00332           return true;
00333         }
00334 
00335       if (rl.lessThan(Expr::handle(l)))
00336         {
00337           if (sign > 0)
00338             {
00339               if (rSign > 0)
00340                 {
00341                   SUNDANCE_VERB_MEDIUM(tab1 << "R is a sum: RL < L < RR, reordering to RL+L+RR");
00342                   rtn = getScalar(rl + Expr::handle(l) + rr);
00343                 }
00344               else
00345                 {
00346                   SUNDANCE_VERB_MEDIUM(tab1 << "R is a sum: RL < L < RR, reordering to RL+L-RR");
00347                   rtn = getScalar(rl + Expr::handle(l) - rr);
00348                 }
00349             }
00350           else
00351             {
00352               if (rSign > 0)
00353                 {
00354                   SUNDANCE_VERB_MEDIUM(tab1 << "R is a sum: RL < L < RR, reordering to -RL+L-RR");
00355                   rtn = getScalar(-rl + Expr::handle(l) - rr);
00356                 }
00357               else
00358                 {
00359                   SUNDANCE_VERB_MEDIUM(tab1 << "R is a sum: RL < L < RR, reordering to -RL+L+RR");
00360                   rtn = getScalar(-rl + Expr::handle(l) + rr);
00361                 }
00362             }
00363           return true;
00364         }
00365     }
00366   
00367   return false;
00368 }
00369 
00370 #endif
00371 
00372 
00373 bool RemoveZeroFromSum::doTransform(const RCP<ScalarExpr>& left, const RCP<ScalarExpr>& right,
00374                                     int sign, RCP<ScalarExpr>& rtn) const
00375 {
00376   TimeMonitor timer(removeZeroTimer());
00377   SUNDANCE_OUT(this->verb() > 1, 
00378                "trying RemoveZeroFromSum");
00379 
00380   /* Check for the trivial case of operation with zero */
00381   
00382   /* If I'm constant and my value is zero, return other */
00383   const ConstantExpr* cl = dynamic_cast<const ConstantExpr*>(left.get());
00384   if (cl != 0 && (cl->value()==0.0 || cl->value()==-0.0))
00385     {
00386       if (verb() > 1)
00387         {
00388           Out::println("RemoveZeroFromSum identified left operand as zero.");
00389           Out::println("Applying transformation 0 + x --> x");
00390         }
00391       rtn = chooseSign(sign, right);
00392       return true;
00393     }
00394 
00395   /* If other is zero, return me */
00396   const ConstantExpr* cr = dynamic_cast<const ConstantExpr*>(right.get());
00397   if (cr != 0 && (cr->value()==0.0 || cr->value()==-0.0)) 
00398     {
00399       if (verb() > 1)
00400         {
00401           Out::println("RemoveZeroFromSum identified right operand as zero.");
00402           Out::println("Applying transformation x + 0 --> x");
00403         }
00404       rtn = left;
00405       return true;
00406     }
00407 
00408   /* otherwise, do no transformation */
00409   return false;
00410 }
00411 
00412 bool MoveConstantsToLeftOfSum::doTransform(const RCP<ScalarExpr>& left, const RCP<ScalarExpr>& right,
00413                                       int sign, RCP<ScalarExpr>& rtn) const
00414 {
00415   TimeMonitor timer(moveConstantsTimer());
00416 
00417   /* if the right operand is a constant, 
00418    * transform u +/- alpha --> +/- alpha + u */
00419   if (right->isConstant())
00420     {
00421       if (verb() > 1)
00422         {
00423           Out::println("MoveConstantsToLeftOfSum identified right "
00424                        "operand as constant.");
00425         }
00426       rtn = getScalar(Expr::handle(chooseSign(sign, right)) 
00427         + Expr::handle(left));
00428       return true;
00429     }
00430 
00431   return false;
00432 }
00433 
00434 
00435 bool RemoveUnaryMinusFromSum::doTransform(const RCP<ScalarExpr>& left,
00436                                           const RCP<ScalarExpr>& right,
00437                                           int sign, 
00438                                           RCP<ScalarExpr>& rtn) const
00439 {
00440   TimeMonitor timer(removeUnaryMinusTimer());
00441   SUNDANCE_OUT(this->verb() > 1, 
00442                "trying RemoveUnaryMinusFromSum");
00443 
00444   /* if the right operand is a unary minus, 
00445    * transform u +/- (-v) --> u -/+ v */
00446   const UnaryMinus* ul = dynamic_cast<const UnaryMinus*>(left.get());
00447   const UnaryMinus* ur = dynamic_cast<const UnaryMinus*>(right.get());
00448   if (ul != 0 && ur != 0)
00449     {
00450       if (verb() > 1)
00451         {
00452           Out::println("RemoveUnaryMinusFromSum identified both "
00453                        "operands as unary minuses.");
00454         }
00455       rtn = getScalar(-(Expr::handle(chooseSign(sign, getScalar(ur->arg()))) 
00456                         + ul->arg()));
00457       return true;
00458     }
00459   else if (ul != 0)
00460     {
00461       if (verb() > 1)
00462         {
00463           Out::println("RemoveUnaryMinusFromSum identified left "
00464                        "operand as unary minus.");
00465         }
00466       if (sign > 0)
00467         {
00468           rtn = getScalar(-(ul->arg() - Expr::handle(right)));
00469         }
00470       else
00471         {
00472           rtn = getScalar(-(ul->arg() + Expr::handle(right)));
00473         }
00474       return true;
00475     }
00476   else if (ur != 0)
00477     {
00478       if (verb() > 1)
00479         {
00480           Out::println("RemoveUnaryMinusFromSum identified right "
00481                        "operand as unary minus.");
00482         }
00483       if (sign > 0)
00484         {
00485           rtn = getScalar(Expr::handle(left) - ur->arg());
00486         }
00487       else
00488         {
00489           rtn = getScalar(Expr::handle(left) + ur->arg());
00490         }
00491       return true;
00492     }
00493 
00494   return false;
00495 }
00496 
00497 bool SumConstants::doTransform(const RCP<ScalarExpr>& left, const RCP<ScalarExpr>& right,
00498                                int sign, RCP<ScalarExpr>& rtn) const
00499 {
00500   
00501   TimeMonitor timer(sumConstantsTimer());
00502   SUNDANCE_OUT(this->verb() > 1, 
00503                "trying SumConstants");
00504 
00505   /* Check to see if both are constant. If so, sum them and return */
00506   if (left->isConstant() && right->isConstant())
00507     {
00508       if (verb() > 1)
00509         {
00510           Out::println("SumConstants identified both "
00511                        "operands as constant. No transformations applied.");
00512         }
00513       rtn = rcp(new SumExpr(left, right, sign));
00514       return true;
00515     }
00516   return false;
00517 }
00518 
00519 bool RearrangeRightSumWithConstant::doTransform(const RCP<ScalarExpr>& left, 
00520                                                 const RCP<ScalarExpr>& right,
00521                                                 int sign, RCP<ScalarExpr>& rtn) const
00522 {
00523   TimeMonitor timer(rearrangeRightSumWithConstantTimer());
00524   const SumExpr* sRight = dynamic_cast<const SumExpr*>(right.get());
00525 
00526   if (sRight != 0)
00527     {
00528       /* The case in which the right operand of right is a constant
00529        * should have been transformed away by now. Do a paranoid 
00530        * check to make sure this hasn't happened */
00531       TEUCHOS_TEST_FOR_EXCEPTION(sRight->rightScalar()->isConstant(),
00532                          std::logic_error,
00533                          "RearrangeRightSumWithConstant: unexpected case, "
00534                          "constant expr"
00535                          << sRight->right() << " found as right operand "
00536                          "in sum " << right->toString());
00537 
00538       if (sRight->leftScalar()->isConstant())
00539         {      
00540           /* If left operand is a constant, transform
00541            * alpha + s1*(beta + s2*u) --> (alpha + s1*beta) + s1*s2*u */
00542           if (left->isConstant())
00543             {
00544               if (verb() > 1)
00545                 {
00546                   Out::println("RearrangeRightSumWithConstant::doTransform "
00547                                "identified right "
00548                                "operand as sum involving a constant, "
00549                                "and left operand as a constant. Applying "
00550                                "transformation alpha + (beta+u) "
00551                                "--> (alpha + beta) + u.");
00552                 }
00553               int s1 = sign;
00554               int s2 = sRight->sign();
00555               Expr alpha = Expr::handle(left);
00556               Expr beta = sRight->left();
00557               Expr u = sRight->right();
00558               rtn = getScalar((alpha + chooseSign(s1, beta)) + chooseSign(s1*s2,u));
00559             }
00560           else  /* if left operand is non-constant, transform
00561                  * u + s1*(alpha + s2*v) --> s1*alpha + (u + s1*s2*v) */
00562             {
00563               if (verb() > 1)
00564                 {
00565                   Out::println("RearrangeRightSumWithConstant::doTransform "
00566                                "identified right "
00567                                "operand as sum involving a constant, "
00568                                "and left operand as non-constant. Applying "
00569                                "transformation u + (alpha + v) "
00570                                "--> alpha + (u + v)");
00571                 }
00572               int s1 = sign;
00573               int s2 = sRight->sign();
00574               Expr u = Expr::handle(left);
00575               Expr alpha = sRight->left();
00576               Expr v = sRight->right();
00577               rtn = getScalar(chooseSign(s1, alpha) + (u + chooseSign(s1*s2, v)));
00578             }
00579           return true;
00580         }
00581     }
00582   return false;
00583 }
00584 
00585 
00586 bool RearrangeLeftSumWithConstant::doTransform(const RCP<ScalarExpr>& left, 
00587                                                 const RCP<ScalarExpr>& right,
00588                                                 int sign, RCP<ScalarExpr>& rtn) const
00589 {
00590   TimeMonitor timer(rearrangeLeftSumWithConstantTimer());
00591   const SumExpr* sLeft = dynamic_cast<const SumExpr*>(left.get());
00592 
00593   if (sLeft != 0 && !left->isConstant())
00594     {
00595       /* The case in which the right operand of left is a constant
00596        * should have been transformed away by now. Do a paranoid 
00597        * check to make sure this hasn't happened */
00598       TEUCHOS_TEST_FOR_EXCEPTION(sLeft->rightScalar()->isConstant(),
00599                          std::logic_error,
00600                          "RearrangeLeftSumWithConstant::doTransform "
00601                          ": unexpected case, constant expr"
00602                          << sLeft->right() << " found as right operand "
00603                          "in sum " << left->toString());
00604       
00605       if (sLeft->leftScalar()->isConstant())
00606         {
00607           /* If right operand is a constant, transform 
00608            * (alpha + s1*u) + s2*beta --> (alpha + s2*beta) + s1*u */
00609           if (right->isConstant())
00610             {
00611               if (verb() > 1)
00612                 {
00613                   Out::println("RearrangeLeftSumWithConstant::doTransform "
00614                                "identified right "
00615                                "operand as constant, "
00616                                "and left operand as sum involving "
00617                                "a constant. Applying "
00618                                "transformation (alpha + u) + beta "
00619                                "--> (alpha + beta) + u.");
00620                 }
00621               int s2 = sign;
00622               int s1 = sLeft->sign();
00623               Expr alpha = sLeft->left();
00624               Expr beta = Expr::handle(right);
00625               Expr u = sLeft->right();
00626               rtn =  getScalar((alpha + chooseSign(s2, beta)) + chooseSign(s1, u));
00627             }
00628           else /* if right operand is a non-constant, transform 
00629                 * (alpha + s1*u) + s2*v --> alpha + (s1*u + s2*v) */
00630             {
00631               if (verb() > 1)
00632                 {
00633                   Out::println("RearrangeLeftSumWithConstant::doTransform "
00634                                "identified right "
00635                                "operand as non-constant, "
00636                                "and left operand as sum involving "
00637                                "a constant. Applying "
00638                                "transformation (alpha + u) + v "
00639                                "--> alpha + (u + v).");
00640                 }
00641               int s2 = sign;
00642               int s1 = sLeft->sign();
00643               Expr alpha = sLeft->left();
00644               Expr u = sLeft->right();
00645               Expr v = Expr::handle(right);
00646               rtn =  getScalar(alpha 
00647                 + (chooseSign(s1, u) + chooseSign(s2, v)));
00648             }
00649           return true;
00650         }
00651     }
00652   return false;
00653 }
00654 
00655 bool SumIntegrals::doTransform(const RCP<ScalarExpr>& left, 
00656                                const RCP<ScalarExpr>& right,
00657                                int sign, RCP<ScalarExpr>& rtn) const
00658 {
00659   TimeMonitor timer(sumIntegralsTimer());
00660   SUNDANCE_OUT(this->verb() > 1, 
00661                "trying SumIntegrals");
00662 
00663   const SumOfIntegrals* sLeft 
00664     = dynamic_cast<const SumOfIntegrals*>(left.get());
00665   const SumOfIntegrals* sRight 
00666     = dynamic_cast<const SumOfIntegrals*>(right.get());
00667 
00668   if (sLeft != 0 || sRight != 0)
00669     {
00670       /* make sure we don't have a case where one is an essential BC and
00671        * the other is not */
00672       bool leftIsBC = (dynamic_cast<const SumOfBCs*>(sLeft) != 0);
00673       bool rightIsBC = (dynamic_cast<const SumOfBCs*>(sRight) != 0);
00674       TEUCHOS_TEST_FOR_EXCEPTION((leftIsBC && !rightIsBC)
00675                          || (!leftIsBC && rightIsBC), std::runtime_error,
00676                          "Attempting to add EssentialBC and non-EssentialBC "
00677                          "integrals: L=" << left->toString() << ", R="
00678                          << right->toString());
00679 
00680       if (sLeft != 0 && sRight != 0)
00681         {
00682           SumOfIntegrals* l;
00683           if (!leftIsBC) l = new SumOfIntegrals(*sLeft);
00684           else l = new SumOfBCs(*dynamic_cast<const SumOfBCs*>(sLeft));
00685           l->merge(sRight, sign);
00686           rtn = rcp(l);
00687           return true;
00688         }
00689 
00690       /* at this point, one of the terms is a global equation. BCs should
00691        * not be involved at this point */
00692       TEUCHOS_TEST_FOR_EXCEPTION(leftIsBC, std::runtime_error,
00693                          "Attempting to add a BC " << left->toString()
00694                          << " and a global expression " << right->toString());
00695 
00696       TEUCHOS_TEST_FOR_EXCEPTION(rightIsBC, std::runtime_error,
00697                          "Attempting to add a BC " << right->toString()
00698                          << " and a global expression " << left->toString());
00699 
00700       if (sLeft != 0 && sRight == 0)
00701         {
00702           SumOfIntegrals* l = new SumOfIntegrals(*sLeft);
00703           const SpatiallyConstantExpr* cRight 
00704             = dynamic_cast<const SpatiallyConstantExpr*>(right.get());
00705 
00706           TEUCHOS_TEST_FOR_EXCEPTION(cRight == 0, std::logic_error,
00707                              "Attempting to add non-constant expression "
00708                              << right->toString() << " to an integral");
00709 
00710           Expr r = Integral(l->nullRegion(), Expr::handle(right));
00711           const SumOfIntegrals* sr 
00712             = dynamic_cast<const SumOfIntegrals*>(r.ptr().get());
00713           l->merge(sr, sign);
00714           rtn = rcp(l);
00715           return true;
00716         }
00717       else
00718         {
00719           SumOfIntegrals* r = new SumOfIntegrals(*sRight);
00720           if (sign < 0) r->changeSign();
00721 
00722           const SpatiallyConstantExpr* cLeft 
00723             = dynamic_cast<const SpatiallyConstantExpr*>(left.get());
00724 
00725           TEUCHOS_TEST_FOR_EXCEPTION(cLeft == 0, std::logic_error,
00726                              "Attempting to add non-constant expression "
00727                              << left->toString() << " to an integral");
00728 
00729           Expr l = Integral(r->nullRegion(), Expr::handle(right));
00730           const SumOfIntegrals* sl 
00731             = dynamic_cast<const SumOfIntegrals*>(l.ptr().get());
00732           r->merge(sl, 1);
00733           rtn = rcp(r);
00734           return true;
00735         }
00736 
00737     }
00738   else
00739     {
00740       return false;
00741     }
00742   TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "this should not happen");
00743   return false; // -Wall;
00744 }
SundanceStdSumTransformations.cpp
