00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceExpr.hpp"
00032 #include "SundanceIntegral.hpp"
00033 #include "SundanceStdProductTransformations.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 #include "SundanceFunctionalDeriv.hpp"
00042 #include "SundanceCoordDeriv.hpp"
00043 #include "SundanceUnaryMinus.hpp"
00044 #include "SundanceZeroExpr.hpp"
00045 #include "SundanceSumOfIntegrals.hpp"
00046 #include "SundanceSymbolicFuncElement.hpp"
00047 #include "SundanceDerivOfSymbFunc.hpp"
00048 #include "SundanceOut.hpp"
00049
00050 using namespace SundanceCore;
00051 using namespace SundanceUtils;
00052
00053 using namespace Teuchos;
00054 using namespace SundanceCore::Internal;
00055
00056
00057 StdProductTransformations::StdProductTransformations()
00058 : ProductTransformationSequence()
00059 {
00060 append(rcp(new RemoveZeroFromProduct()));
00061 append(rcp(new RemoveOneFromProduct()));
00062 append(rcp(new RemoveMinusOneFromProduct()));
00063 append(rcp(new KillDiffOpOnConstant()));
00064 append(rcp(new BringConstantOutsideDiffOp()));
00065 append(rcp(new MoveUnaryMinusOutsideProduct()));
00066 append(rcp(new AssociateHungryDiffOpWithOperand()));
00067 append(rcp(new DistributeSumOfDiffOps()));
00068 append(rcp(new ApplySimpleDiffOp()));
00069 append(rcp(new TakeConstantUnderIntegralSign()));
00070 append(rcp(new MultiplyConstants()));
00071 append(rcp(new MoveConstantsToLeftOfProduct()));
00072 append(rcp(new RearrangeRightProductWithConstant()));
00073 append(rcp(new RearrangeLeftProductWithConstant()));
00074 }
00075
00076 bool RemoveZeroFromProduct::doTransform(const RefCountPtr<ScalarExpr>& left,
00077 const RefCountPtr<ScalarExpr>& right,
00078 RefCountPtr<ScalarExpr>& rtn) const
00079 {
00080 SUNDANCE_OUT(this->verbosity() > VerbLow,
00081 "trying RemoveZerofromProduct");
00082
00083
00084 const ConstantExpr* cl = dynamic_cast<const ConstantExpr*>(left.get());
00085 const ConstantExpr* cr = dynamic_cast<const ConstantExpr*>(right.get());
00086
00087 if (cl != 0)
00088 {
00089 if (cl->value()==0.0 || cl->value()==-0.0)
00090 {
00091 if (verbosity() > 1)
00092 {
00093 Out::println("RemoveOneFromProduct::doTransform "
00094 "identified multiplication "
00095 "by zero. Applying transformation 0*u --> 0");
00096 }
00097 rtn = rcp(new ZeroExpr());
00098 return true;
00099 }
00100 }
00101 if (cr != 0)
00102 {
00103 if (cr->value()==0.0 || cr->value()==-0.0)
00104 {
00105 if (verbosity() > 1)
00106 {
00107 Out::println("RemoveOneFromProduct::doTransform "
00108 "identified multiplication "
00109 "by zero. Applying transformation u*0 --> u");
00110 }
00111 rtn = rcp(new ZeroExpr());
00112 return true;
00113 }
00114 }
00115 return false;
00116 }
00117
00118 bool RemoveOneFromProduct::doTransform(const RefCountPtr<ScalarExpr>& left,
00119 const RefCountPtr<ScalarExpr>& right,
00120 RefCountPtr<ScalarExpr>& rtn) const
00121 {
00122 SUNDANCE_OUT(this->verbosity() > VerbLow,
00123 "trying RemoveOnefromProduct");
00124
00125 const ConstantExpr* cl = dynamic_cast<const ConstantExpr*>(left.get());
00126 const ConstantExpr* cr = dynamic_cast<const ConstantExpr*>(right.get());
00127
00128 if (cl != 0)
00129 {
00130 if (cl->value()==1.0)
00131 {
00132 if (verbosity() > 1)
00133 {
00134 Out::println("RemoveOneFromProduct::doTransform "
00135 "identified multiplication "
00136 "by one. Applying transformation 1*u --> u");
00137 }
00138 rtn = getScalar(Expr::handle(right));
00139 return true;
00140 }
00141 }
00142 if (cr != 0)
00143 {
00144 if (cr->value()==1.0)
00145 {
00146 if (verbosity() > 1)
00147 {
00148 Out::println("RemoveOneFromProduct::doTransform "
00149 "identified multiplication "
00150 "by one. Applying transformation u*1 --> u");
00151 }
00152 rtn = getScalar(Expr::handle(left));
00153 return true;
00154 }
00155 }
00156 return false;
00157 }
00158
00159
00160 bool RemoveMinusOneFromProduct::doTransform(const RefCountPtr<ScalarExpr>& left,
00161 const RefCountPtr<ScalarExpr>& right,
00162 RefCountPtr<ScalarExpr>& rtn) const
00163 {
00164 SUNDANCE_OUT(this->verbosity() > VerbLow,
00165 "trying RemoveOnefromProduct");
00166
00167 const ConstantExpr* cl = dynamic_cast<const ConstantExpr*>(left.get());
00168 const ConstantExpr* cr = dynamic_cast<const ConstantExpr*>(right.get());
00169
00170 if (cl != 0)
00171 {
00172 if (cl->value()==-1.0)
00173 {
00174 if (verbosity() > 1)
00175 {
00176 Out::println("RemoveMinusOneFromProduct::doTransform "
00177 "identified multiplication "
00178 "by one. Applying transformation -1*u --> -u");
00179 }
00180 rtn = getScalar(-Expr::handle(right));
00181 return true;
00182 }
00183 }
00184 if (cr != 0)
00185 {
00186 if (cr->value()==-1.0)
00187 {
00188 if (verbosity() > 1)
00189 {
00190 Out::println("RemoveMinusOneFromProduct::doTransform "
00191 "identified multiplication "
00192 "by one. Applying transformation u*(-1) --> -u");
00193 }
00194 rtn = getScalar(-Expr::handle(left));
00195 return true;
00196 }
00197 }
00198 return false;
00199 }
00200
00201 bool MoveConstantsToLeftOfProduct::doTransform(const RefCountPtr<ScalarExpr>& left,
00202 const RefCountPtr<ScalarExpr>& right,
00203 RefCountPtr<ScalarExpr>& rtn) const
00204 {
00205 SUNDANCE_OUT(this->verbosity() > VerbLow,
00206 "trying MoveConstantsToLeftOfProduct");
00207
00208
00209
00210 if (!left->isConstant() && right->isConstant())
00211 {
00212 if (verbosity() > 1)
00213 {
00214 Out::println("MoveConstantsToLeftOfProduct::doTransform "
00215 "identified right operand "
00216 "as a constant. Applying transformation u*alpha "
00217 "--> alpha*u.");
00218 }
00219 rtn = getScalar(Expr::handle(right) * Expr::handle(left));
00220 return true;
00221 }
00222 return false;
00223 }
00224
00225 bool MoveUnaryMinusOutsideProduct::doTransform(const RefCountPtr<ScalarExpr>& left,
00226 const RefCountPtr<ScalarExpr>& right,
00227 RefCountPtr<ScalarExpr>& rtn) const
00228 {
00229 SUNDANCE_OUT(this->verbosity() > VerbLow,
00230 "trying MoveUnaryMinusOutsideProduct");
00231
00232
00233
00234 const UnaryMinus* ul = dynamic_cast<const UnaryMinus*>(left.get());
00235 const UnaryMinus* ur = dynamic_cast<const UnaryMinus*>(right.get());
00236 if (ur != 0 && ul != 0)
00237 {
00238 if (verbosity() > 1)
00239 {
00240 Out::println("MoveUnaryMinusOutsideProduct::doTransform "
00241 "identified both operands "
00242 "as unary minuses. Applying transformation (-x)*(-y) "
00243 "--> x*y.");
00244 }
00245 rtn = getScalar(ul->arg() * ur->arg());
00246 return true;
00247 }
00248 else if (ur != 0)
00249 {
00250 if (verbosity() > 1)
00251 {
00252 Out::println("MoveUnaryMinusOutsideProduct::doTransform "
00253 "identified right operand "
00254 "as a unary minus. Applying transformation x*(-y) "
00255 "--> -(x*y).");
00256 }
00257 rtn = rcp(new UnaryMinus(getScalar(Expr::handle(left) * ur->arg())));
00258 return true;
00259 }
00260 else if (ul != 0)
00261 {
00262 if (verbosity() > 1)
00263 {
00264 Out::println("MoveUnaryMinusOutsideProduct::doTransform "
00265 "identified left operand "
00266 "as a unary minus. Applying transformation (-x)*y "
00267 "--> -(x*y).");
00268 }
00269 rtn = rcp(new UnaryMinus(getScalar(ul->arg() *Expr::handle(right))));
00270 return true;
00271 }
00272 return false;
00273 }
00274
00275 bool MultiplyConstants::doTransform(const RefCountPtr<ScalarExpr>& left, const RefCountPtr<ScalarExpr>& right,
00276 RefCountPtr<ScalarExpr>& rtn) const
00277 {
00278 SUNDANCE_OUT(this->verbosity() > VerbLow,
00279 "trying MultiplyConstants");
00280
00281 if (left->isConstant() && right->isConstant())
00282 {
00283 if (left->isImmutable() && right->isImmutable())
00284 {
00285 const ConstantExpr* cl = dynamic_cast<const ConstantExpr*>(left.get());
00286 const ConstantExpr* cr = dynamic_cast<const ConstantExpr*>(right.get());
00287 TEST_FOR_EXCEPTION(cl==0 || cr==0, InternalError,
00288 "MultiplyConstants::doTransform() logic error: "
00289 "L and R identified as immutable, but could "
00290 "not be cast to ConstantExprs");
00291 rtn = rcp(new ConstantExpr(cl->value() * cr->value()));
00292 return true;
00293 }
00294 if (verbosity() > 1)
00295 {
00296 Out::println("MultiplyConstants::doTransform "
00297 "identified both operands "
00298 "as constants. Forming the product without any "
00299 "transformation");
00300 }
00301 rtn = rcp(new ProductExpr(left, right));
00302 return true;
00303 }
00304 return false;
00305 }
00306
00307 bool AssociateHungryDiffOpWithOperand::doTransform(const RefCountPtr<ScalarExpr>& left,
00308 const RefCountPtr<ScalarExpr>& right,
00309 RefCountPtr<ScalarExpr>& rtn) const
00310 {
00311 SUNDANCE_OUT(this->verbosity() > VerbLow,
00312 "trying AssociateHungryDiffOpWithOperand");
00313 if (left->isHungryDiffOp())
00314 {
00315
00316
00317
00318 const ProductExpr* pLeft
00319 = dynamic_cast<const ProductExpr*>(left.get());
00320 if (pLeft != 0)
00321 {
00322 Expr ll = pLeft->left();
00323 Expr lr = pLeft->right();
00324 if (verbosity() > 1)
00325 {
00326 Out::println("AssociateHungryDiffOpWithOperand::doTransform "
00327 "identified left "
00328 "operand as a product with a hungry diff op "
00329 "as the last factor. "
00330 "Applying (u*D)*v --> u*(D*v).");
00331 }
00332 rtn = getScalar(ll*(lr*Expr::handle(right)));
00333 return true;
00334 }
00335 }
00336 return false;
00337 }
00338
00339 bool KillDiffOpOnConstant::doTransform(const RefCountPtr<ScalarExpr>& left,
00340 const RefCountPtr<ScalarExpr>& right,
00341 RefCountPtr<ScalarExpr>& rtn) const
00342 {
00343 SUNDANCE_OUT(this->verbosity() > VerbLow, "trying KillDiffOpOnConstant");
00344 if (left->isHungryDiffOp())
00345 {
00346
00347
00348 if (right->isConstant())
00349 {
00350 rtn = rcp(new ZeroExpr());
00351 if (verbosity() > 1)
00352 {
00353 Out::println("KillDiffOpOnConstant::doTransform "
00354 "identified constant "
00355 "as operand of diff op. Applying D*alpha --> 0");
00356 }
00357 return true;
00358 }
00359
00360
00361 const SumExpr* sRight = dynamic_cast<const SumExpr*>(right.get());
00362 if (sRight != 0 && sRight->leftScalar()->isConstant())
00363 {
00364 if (verbosity() > 1)
00365 {
00366 Out::println("KillDiffOpOnConstant::doTransform "
00367 "identified constant "
00368 "term in operand of diff op. "
00369 "Applying D*(alpha+u) --> D*u");
00370 }
00371 rtn = getScalar(chooseSign(sRight->sign(), Expr::handle(left)) * sRight->right());
00372 return true;
00373 }
00374 }
00375 return false;
00376 }
00377
00378 bool BringConstantOutsideDiffOp::doTransform(const RefCountPtr<ScalarExpr>& left,
00379 const RefCountPtr<ScalarExpr>& right,
00380 RefCountPtr<ScalarExpr>& rtn) const
00381 {
00382
00383 SUNDANCE_OUT(this->verbosity() > VerbLow, "trying BringConstantOutsideDiffOp");
00384 if (left->isHungryDiffOp())
00385 {
00386
00387 const ProductExpr* pRight
00388 = dynamic_cast<const ProductExpr*>(right.get());
00389 if (pRight != 0 && pRight->leftScalar()->isConstant())
00390 {
00391 if (verbosity() > 1)
00392 {
00393 Out::println("BringConstantOutsideDiffOp::doTransform "
00394 "identified constant "
00395 "coefficient in operand of diff op. "
00396 "Applying D*(alpha*u) --> alpha*D*u");
00397 }
00398 rtn = getScalar(pRight->left() * (Expr::handle(left) * pRight->right()));
00399 return true;
00400 }
00401 }
00402 return false;
00403 }
00404
00405 bool DistributeSumOfDiffOps::doTransform(const RefCountPtr<ScalarExpr>& left,
00406 const RefCountPtr<ScalarExpr>& right,
00407 RefCountPtr<ScalarExpr>& rtn) const
00408 {
00409 SUNDANCE_OUT(this->verbosity() > VerbLow,"trying DistributeSumOfDiffOps");
00410 if (left->isHungryDiffOp())
00411 {
00412
00413
00414 const SumExpr* sLeft = dynamic_cast<const SumExpr*>(left.get());
00415 if (sLeft != 0)
00416 {
00417 Expr ll = sLeft->left();
00418 Expr lr = sLeft->right();
00419 if (verbosity() > 1)
00420 {
00421 Out::println("DistributeSumOfDiffOps::doTransform "
00422 "identified left "
00423 "operand as a sum of hungry diff ops. "
00424 "Applying (D1 + D2)*u --> D1*u + D2*u");
00425 }
00426 rtn = getScalar(ll*Expr::handle(right) + chooseSign(sLeft->sign(), lr)*Expr::handle(right));
00427 return true;
00428 }
00429 }
00430 return false;
00431 }
00432
00433 bool ApplySimpleDiffOp::doTransform(const RefCountPtr<ScalarExpr>& left,
00434 const RefCountPtr<ScalarExpr>& right,
00435 RefCountPtr<ScalarExpr>& rtn) const
00436 {
00437 SUNDANCE_OUT(this->verbosity() > VerbLow, "trying ApplySimpleDiffOp");
00438 if (left->isHungryDiffOp())
00439 {
00440 const Derivative* dLeft
00441 = dynamic_cast<const Derivative*>(left.get());
00442 if (dLeft != 0)
00443 {
00444 const SymbolicFuncElement* sf
00445 = dynamic_cast<const SymbolicFuncElement*>(right.get());
00446 if (sf != 0 && optimizeFunctionDiffOps())
00447 {
00448 rtn = rcp(new DerivOfSymbFunc(dLeft->multiIndex(), right));
00449 }
00450 else
00451 {
00452 rtn = rcp(new DiffOp(dLeft->multiIndex(), right));
00453 }
00454 return true;
00455 }
00456 }
00457 return false;
00458 }
00459
00460 bool RearrangeRightProductWithConstant::doTransform(const RefCountPtr<ScalarExpr>& left,
00461 const RefCountPtr<ScalarExpr>& right,
00462 RefCountPtr<ScalarExpr>& rtn) const
00463 {
00464
00465
00466 const ProductExpr* pRight
00467 = dynamic_cast<const ProductExpr*>(right.get());
00468
00469
00470
00471
00472
00473 TEST_FOR_EXCEPTION(pRight != 0 && pRight->rightScalar()->isConstant(),
00474 InternalError,
00475 "unexpected case in "
00476 "RearrangeRightProductWithConstant::doTransform: "
00477 "the right operand "
00478 << pRight->right()
00479 << "of the right operand " << right->toString()
00480 << " is a constant. That case should have been "
00481 "transformed out by now.");
00482
00483 if (pRight != 0 && !pRight->isConstant()
00484 && pRight->leftScalar()->isConstant())
00485 {
00486
00487
00488
00489 if (left->isConstant())
00490 {
00491 if (verbosity() > 1)
00492 {
00493 Out::println("RearrangeRightProductWithConstant::doTransform: "
00494 "identified left operand "
00495 "as a constant, and right operand as a product "
00496 "involving a constant. Applying transformation "
00497 "alpha*(beta*u) --> (alpha*beta)*u");
00498 }
00499 rtn = getScalar((Expr::handle(left) * pRight->left()) * pRight->right());
00500 return true;
00501 }
00502 else
00503
00504
00505
00506 {
00507 if (verbosity() > 1)
00508 {
00509 Out::println("RearrangeRightProductWithConstant::doTransform: "
00510 "identified left operand "
00511 "as non-constant, and right operand as a product "
00512 "involving a constant. Applying transformation "
00513 "u * (alpha*v) --> alpha*(u*v)");
00514 }
00515 rtn = getScalar(pRight->left() * (Expr::handle(left) * pRight->right()));
00516 return true;
00517 }
00518 }
00519 return false;
00520 }
00521
00522
00523 bool RearrangeLeftProductWithConstant::doTransform(const RefCountPtr<ScalarExpr>& left,
00524 const RefCountPtr<ScalarExpr>& right,
00525 RefCountPtr<ScalarExpr>& rtn) const
00526 {
00527
00528
00529
00530
00531
00532 const ProductExpr* pLeft
00533 = dynamic_cast<const ProductExpr*>(left.get());
00534
00535 if (pLeft != 0 && !pLeft->isConstant()
00536 && pLeft->leftScalar()->isConstant())
00537 {
00538
00539
00540
00541 TEST_FOR_EXCEPTION(pLeft != 0 && pLeft->rightScalar()->isConstant(),
00542 InternalError,
00543 "RearrangeLeftProductWithConstant::doTransform: "
00544 "the right operand "
00545 << pLeft->right()
00546 << "of the left operand " << left->toString()
00547 << " is a constant. That case should have been "
00548 "transformed out by now.");
00549
00550
00551 if (right->isConstant())
00552 {
00553 if (verbosity() > 1)
00554 {
00555 Out::println("RearrangeLeftProductWithConstant::doTransform: "
00556 "identified right operand "
00557 "as a constant, and left operand as a product "
00558 "involving a constant. Applying transformation "
00559 "(alpha*u)*beta --> (alpha*beta)*u");
00560 }
00561 rtn = getScalar((pLeft->left() * Expr::handle(right)) * pLeft->right());
00562 return true;
00563 }
00564 else
00565
00566
00567 {
00568 if (verbosity() > 1)
00569 {
00570 Out::println("RearrangeLeftProductWithConstant::doTransform: "
00571 "identified right operand "
00572 "as non-constant, and left operand as a product "
00573 "involving a constant. Applying transformation "
00574 "(alpha*u)*v --> alpha*(u*v)");
00575 }
00576 rtn = getScalar(pLeft->left() * (pLeft->right() * Expr::handle(right)));
00577 return true;
00578 }
00579 }
00580 return false;
00581 }
00582
00583
00584 bool TakeConstantUnderIntegralSign::doTransform(const RefCountPtr<ScalarExpr>& left,
00585 const RefCountPtr<ScalarExpr>& right,
00586 RefCountPtr<ScalarExpr>& rtn) const
00587 {
00588 const SumOfIntegrals* sLeft
00589 = dynamic_cast<const SumOfIntegrals*>(left.get());
00590 const SumOfIntegrals* sRight
00591 = dynamic_cast<const SumOfIntegrals*>(right.get());
00592
00593 TEST_FOR_EXCEPTION(sLeft != 0 && sRight != 0, InternalError,
00594 "Product of integrals detected: left="
00595 << left->toString() << " right=" << right->toString());
00596
00597 if (sLeft != 0 || sRight != 0)
00598 {
00599 if (sLeft != 0)
00600 {
00601 SumOfIntegrals* l = new SumOfIntegrals(*sLeft);
00602 const SpatiallyConstantExpr* cRight
00603 = dynamic_cast<const SpatiallyConstantExpr*>(right.get());
00604 TEST_FOR_EXCEPTION(cRight == 0, InternalError,
00605 "Attempting to multiply non-constant expression "
00606 << right->toString() << " with an integral");
00607 l->multiplyByConstant(cRight);
00608
00609 rtn = rcp(l);
00610 return true;
00611 }
00612
00613 if (sRight != 0)
00614 {
00615 SumOfIntegrals* r = new SumOfIntegrals(*sRight);
00616 const SpatiallyConstantExpr* cLeft
00617 = dynamic_cast<const SpatiallyConstantExpr*>(left.get());
00618 TEST_FOR_EXCEPTION(cLeft == 0, InternalError,
00619 "Attempting to multiply non-constant expression "
00620 << left->toString() << " with an integral");
00621 r->multiplyByConstant(cLeft);
00622
00623 rtn = rcp(r);
00624 return true;
00625 }
00626 }
00627 else
00628 {
00629 return false;
00630 }
00631 return false;
00632 }