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