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 "SundanceListExpr.hpp"
00033 #include "SundanceSumExpr.hpp"
00034 #include "SundanceProductExpr.hpp"
00035 #include "SundanceConstantExpr.hpp"
00036 #include "SundanceCoordExpr.hpp"
00037 #include "SundanceDerivative.hpp"
00038 #include "SundanceDiffOp.hpp"
00039
00040
00041 #include "SundanceUnaryMinus.hpp"
00042 #include "SundanceZeroExpr.hpp"
00043 #include "SundanceComplexExpr.hpp"
00044 #include "SundanceOut.hpp"
00045 #include "PlayaTabs.hpp"
00046 #include "SundanceStdSumTransformations.hpp"
00047 #include "SundanceStdProductTransformations.hpp"
00048 #include "SundanceNonlinearUnaryOp.hpp"
00049 #include "SundanceStdMathOps.hpp"
00050 #include "SundanceParameter.hpp"
00051 #include "SundanceSpectralBasis.hpp"
00052 #include "SundanceSpectralExpr.hpp"
00053 #include "SundanceUnknownFuncElement.hpp"
00054 #include "SundanceTestFuncElement.hpp"
00055 #include "Teuchos_TimeMonitor.hpp"
00056
00057 using namespace Sundance;
00058 using namespace Sundance;
00059
00060 using namespace Teuchos;
00061 using namespace Sundance;
00062
00063 static Time& sumTimer()
00064 {
00065 static RCP<Time> rtn
00066 = TimeMonitor::getNewTimer("symbolic sum");
00067 return *rtn;
00068 }
00069
00070 static Time& unaryMinusTimer()
00071 {
00072 static RCP<Time> rtn
00073 = TimeMonitor::getNewTimer("symbolic unary minus");
00074 return *rtn;
00075 }
00076
00077
00078 static Time& trySumComplexTimer()
00079 {
00080 static RCP<Time> rtn
00081 = TimeMonitor::getNewTimer("test for complex sum");
00082 return *rtn;
00083 }
00084
00085 static Time& tryMultiplyComplexTimer()
00086 {
00087 static RCP<Time> rtn
00088 = TimeMonitor::getNewTimer("test for complex product");
00089 return *rtn;
00090 }
00091
00092
00093
00094 static Time& prodTimer()
00095 {
00096 static RCP<Time> rtn
00097 = TimeMonitor::getNewTimer("symbolic product");
00098 return *rtn;
00099 }
00100
00101
00102 Expr::Expr(const double& c)
00103 : Playa::Handle<ExprBase>(new ConstantExpr(c))
00104 {}
00105
00106 Expr::Expr(const std::complex<double>& c)
00107 : Playa::Handle<ExprBase>(new ComplexExpr(new ConstantExpr(c.real()),
00108 new ConstantExpr(c.imag())))
00109 {}
00110
00111 bool Expr::isComplex() const
00112 {
00113 return dynamic_cast<const ComplexExpr*>(ptr().get()) != 0;
00114 }
00115
00116
00117
00118 bool Expr::isSpectral() const
00119 {
00120 return dynamic_cast<const SpectralExpr*>(ptr().get()) != 0;
00121 }
00122
00123
00124
00125
00126 XMLObject Expr::toXML() const
00127 {
00128 TimeMonitor t(outputTimer());
00129
00130 return ptr()->toXML();
00131 }
00132
00133 string Expr::toString() const
00134 {
00135 TimeMonitor t(outputTimer());
00136
00137 TeuchosOStringStream ss;
00138 ptr()->toText(ss, false);
00139 return TEUCHOS_OSTRINGSTREAM_GET_C_STR(ss);
00140 }
00141
00142
00143 bool Expr::sameAs(const Expr& other) const
00144 {
00145 if (this->lessThan(other)) return false;
00146 if (other.lessThan(*this)) return false;
00147 return true;
00148 }
00149
00150 bool Expr::operator<(const Expr& other) const
00151 {
00152 return this->lessThan(other);
00153 }
00154
00155 bool Expr::lessThan(const Expr& other) const
00156 {
00157 RCP<ScalarExpr> sThis = rcp_dynamic_cast<ScalarExpr>(ptr());
00158 RCP<ScalarExpr> sOther = rcp_dynamic_cast<ScalarExpr>(other.ptr());
00159
00160 TEUCHOS_TEST_FOR_EXCEPTION(sThis.get()==0, std::logic_error,
00161 "ordering not defined for non-scalar expression "
00162 << toString());
00163
00164 TEUCHOS_TEST_FOR_EXCEPTION(sOther.get()==0, std::logic_error,
00165 "ordering not defined for non-scalar expressions"
00166 << other.toString());
00167
00168
00169 const ConstantExpr* cMe = dynamic_cast<const ConstantExpr*>(sThis.get());
00170 const ConstantExpr* cOther = dynamic_cast<const ConstantExpr*>(sOther.get());
00171
00172
00173 if (cMe != 0 && cOther==0) return true;
00174 if (cOther != 0 && cMe==0) return false;
00175 if (cOther != 0 && cMe != 0) return cMe->lessThan(cOther);
00176
00177
00178
00179 const SpatiallyConstantExpr* scMe
00180 = dynamic_cast<const SpatiallyConstantExpr*>(sThis.get());
00181 const SpatiallyConstantExpr* scOther
00182 = dynamic_cast<const SpatiallyConstantExpr*>(sOther.get());
00183 if (scMe != 0 && scOther==0) return true;
00184 if (scOther != 0 && scMe==0) return false;
00185
00186
00187
00188 if (sThis->typeName() < sOther->typeName()) return true;
00189 if (sThis->typeName() > sOther->typeName()) return false;
00190
00191
00192 return sThis->lessThan(sOther.get());
00193 }
00194
00195 Sundance::Map<Expr, int> Expr::getSumTree() const
00196 {
00197 RCP<ScalarExpr> sThis = rcp_dynamic_cast<ScalarExpr>(ptr());
00198
00199 TEUCHOS_TEST_FOR_EXCEPTION(sThis.get()==0, std::logic_error,
00200 "etSumTree() not defined for non-scalar expression "
00201 << toString());
00202
00203 const SumExpr* s = dynamic_cast<const SumExpr*>(sThis.get());
00204 const UnaryMinus* u = dynamic_cast<const UnaryMinus*>(sThis.get());
00205 if (s != 0)
00206 {
00207 return s->getSumTree();
00208 }
00209 else if (u != 0)
00210 {
00211 Sundance::Map<Expr, int> rtn;
00212 rtn.put(u->arg(), -1);
00213 return rtn;
00214 }
00215 else
00216 {
00217 Sundance::Map<Expr, int> rtn;
00218 rtn.put(*this, 1);
00219 return rtn;
00220 }
00221
00222 }
00223
00224 bool Expr::tryAddComplex(const Expr& L, const Expr& R, int sign,
00225 Expr& rtn) const
00226 {
00227 TimeMonitor t(trySumComplexTimer());
00228 TEUCHOS_TEST_FOR_EXCEPTION(L.size() != 1 || R.size() != 1, std::logic_error,
00229 "non-scalar exprs should have been reduced before "
00230 "call to tryAddComplex(). Left=" << L << " right="
00231 << R);
00232 if (L.isComplex() || R.isComplex())
00233 {
00234 if (sign > 0)
00235 {
00236 rtn = new ComplexExpr(L.real() + R.real(),
00237 L.imag() + R.imag());
00238 }
00239 else
00240 {
00241 rtn = new ComplexExpr(L.real() - R.real(),
00242 L.imag() - R.imag());
00243 }
00244 return true;
00245 }
00246 else
00247 {
00248 return false;
00249 }
00250 }
00251
00252
00253
00254 bool Expr::tryMultiplyComplex(const Expr& L, const Expr& R,
00255 Expr& rtn) const
00256 {
00257 TimeMonitor t(tryMultiplyComplexTimer());
00258 TEUCHOS_TEST_FOR_EXCEPTION(L.size() != 1 || R.size() != 1, std::logic_error,
00259 "non-scalar exprs should have been reduced before "
00260 "call to tryMultiplyComplex(). Left=" << L << " right="
00261 << R);
00262
00263 if (L.isComplex() || R.isComplex())
00264 {
00265 if (Re(L).sameAs(Re(R)) && Im(L).sameAs(-Im(R)))
00266 {
00267 rtn = Re(R)*Re(R) + Im(R)*Im(R);
00268 }
00269 else
00270 {
00271 Expr re = Re(L)*Re(R) - Im(L)*Im(R);
00272 Expr im = Re(L)*Im(R) + Im(L)*Re(R);
00273 rtn = new ComplexExpr(re, im);
00274 }
00275 return true;
00276 }
00277 else
00278 {
00279 return false;
00280 }
00281 }
00282
00283
00284
00285
00286 Expr Expr::operator+(const Expr& other) const
00287 {
00288 TimeMonitor t(opTimer());
00289 TimeMonitor ts(sumTimer());
00290
00291
00292 if (this->size()!=1 || other.size()!=1)
00293 {
00294 TEUCHOS_FUNC_TIME_MONITOR("plus branch 1");
00295 Array<Expr> rtn(this->size());
00296 TEUCHOS_TEST_FOR_EXCEPTION(this->size() != other.size(), std::runtime_error,
00297 "mismatched list structures in operands L="
00298 << *this << ", R=" << other);
00299 for (int i=0; i<this->size(); i++)
00300 {
00301 rtn[i] = (*this)[i] + other[i];
00302 }
00303 return new ListExpr(rtn);
00304 }
00305 else
00306 {
00307 TEUCHOS_FUNC_TIME_MONITOR("plus branch 2");
00308 Expr rtn;
00309
00310 if (tryAddComplex((*this)[0], other[0], 1, rtn)) return rtn;
00311
00312 return (*this)[0].sum(other[0], 1);
00313 }
00314 }
00315
00316 Expr Expr::operator-(const Expr& other) const
00317 {
00318 TimeMonitor t(opTimer());
00319 TimeMonitor ts(sumTimer());
00320
00321
00322 if (this->size()!=1 || other.size()!=1)
00323 {
00324 TEUCHOS_FUNC_TIME_MONITOR("minus branch 1");
00325 Array<Expr> rtn(this->size());
00326 TEUCHOS_TEST_FOR_EXCEPTION(this->size() != other.size(), std::runtime_error,
00327 "mismatched list structures in operands L="
00328 << *this << ", R=" << other);
00329 for (int i=0; i<this->size(); i++)
00330 {
00331 rtn[i] = (*this)[i] - other[i];
00332 }
00333 return new ListExpr(rtn);
00334 }
00335 else
00336 {
00337 TEUCHOS_FUNC_TIME_MONITOR("minus branch 2");
00338 Expr rtn;
00339
00340 if (tryAddComplex((*this)[0], other[0], -1, rtn)) return rtn;
00341
00342 return (*this)[0].sum(other[0], -1);
00343 }
00344 }
00345
00346
00347 Expr Expr::sum(const Expr& other, int sign) const
00348 {
00349 TEUCHOS_FUNC_TIME_MONITOR("Expr::sum()");
00350 RCP<ScalarExpr> rtn;
00351 RCP<ScalarExpr> sThis = rcp_dynamic_cast<ScalarExpr>(ptr());
00352 RCP<ScalarExpr> sOther = rcp_dynamic_cast<ScalarExpr>(other.ptr());
00353
00354 TEUCHOS_TEST_FOR_EXCEPTION(ptr().get()==NULL, std::logic_error,
00355 "Expr::sum() detected null this pointer");
00356
00357 TEUCHOS_TEST_FOR_EXCEPTION(sThis.get()==NULL, std::logic_error,
00358 "Expr::sum(): Left operand " << toString()
00359 << " is a non-scalar expression. All list structure "
00360 "should have been handled before this point");
00361
00362 TEUCHOS_TEST_FOR_EXCEPTION(sOther.get()==NULL, std::logic_error,
00363 "Expr::sum(): Right operand " << other.toString()
00364 << " is a non-scalar expression. All list structure "
00365 "should have been handled before this point");
00366
00367 static StdSumTransformations trans;
00368
00369 if (trans.doTransform(sThis, sOther, sign, rtn))
00370 {
00371 if (SymbolicTransformation::classVerbosity() > 0)
00372 {
00373 Out::println("Expr::sum() transformed sum\n["
00374 + toString() + "+"
00375 + other.toString() + "]\n to\n ["
00376 + rtn->toString() + "]");
00377 }
00378 return handle(rtn);
00379 }
00380
00381 else return new SumExpr(sThis, sOther, sign);
00382 }
00383
00384
00385 Expr Expr::operator*(const Expr& other) const
00386 {
00387 TimeMonitor t(opTimer());
00388 TimeMonitor tp(prodTimer());
00389 Tabs tab1;
00390
00391
00392 if (this->size() == 1 && other.size()==1)
00393 {
00394 Expr rtn;
00395
00396 if (tryMultiplyComplex((*this)[0], other[0], rtn))
00397 {
00398 return rtn;
00399 }
00400
00401 rtn = (*this)[0].multiply(other[0]);
00402 return rtn;
00403 }
00404
00405
00406 if (this->size()==1)
00407 {
00408 Array<Expr> rtn(other.size());
00409 for (int i=0; i<other.size(); i++)
00410 {
00411 rtn[i] = (*this)[0] * other[i];
00412 }
00413 return new ListExpr(rtn);
00414 }
00415
00416
00417 if (other.size()==1)
00418 {
00419 Array<Expr> rtn(this->size());
00420 for (int i=0; i<this->size(); i++)
00421 {
00422 rtn[i] = (*this)[i] * other[0];
00423 }
00424 return new ListExpr(rtn);
00425 }
00426
00427
00428 if (this->size() == totalSize() && other.size()==other.totalSize()
00429 && this->size() == other.size())
00430 {
00431 Expr rtn = new ZeroExpr();
00432
00433 for (int i=0; i<this->size(); i++)
00434 {
00435 rtn = rtn + (*this)[i]*other[i];
00436 }
00437 return rtn;
00438 }
00439
00440
00441
00442 int cols = (*this)[0].size();
00443 bool rectangular = true;
00444 for (int i=0; i<this->size(); i++)
00445 {
00446 if ((*this)[i].size() != cols) rectangular = false;
00447 }
00448 TEUCHOS_TEST_FOR_EXCEPTION(!rectangular, std::runtime_error,
00449 "Expr::operator* detected list-list multiplication "
00450 "with a non-rectangular left operator "
00451 << toString());
00452
00453 TEUCHOS_TEST_FOR_EXCEPTION(cols != other.size(), std::runtime_error,
00454 "Expr::operator* detected mismatched dimensions in "
00455 "list-list multiplication. Left operator is "
00456 << toString() << ", right operator is "
00457 << other.toString());
00458
00459 Array<Expr> rtn(this->size());
00460 for (int i=0; i<this->size(); i++)
00461 {
00462 rtn[i] = (*this)[i] * other;
00463 }
00464
00465 return new ListExpr(rtn);
00466 }
00467
00468 Expr Expr::divide(const Expr& other) const
00469 {
00470 RCP<ScalarExpr> sOther = rcp_dynamic_cast<ScalarExpr>(other.ptr());
00471 Expr recip = new NonlinearUnaryOp(sOther, rcp(new StdReciprocal()));
00472 return (*this)[0].multiply(recip);
00473 }
00474
00475 Expr Expr::multiply(const Expr& other) const
00476 {
00477 RCP<ScalarExpr> rtn;
00478 RCP<ScalarExpr> sThis = rcp_dynamic_cast<ScalarExpr>(ptr());
00479 RCP<ScalarExpr> sOther = rcp_dynamic_cast<ScalarExpr>(other.ptr());
00480
00481 TEUCHOS_TEST_FOR_EXCEPTION(sThis.get()==NULL, std::logic_error,
00482 "Expr::multiply(): Left operand " << toString()
00483 << " is a non-scalar expression. All list structure "
00484 "should have been handled before this point");
00485
00486 TEUCHOS_TEST_FOR_EXCEPTION(sOther.get()==NULL, std::logic_error,
00487 "Expr::multiply(): Right operand " << other.toString()
00488 << " is a non-scalar expression. All list structure "
00489 "should have been handled before this point");
00490
00491 static StdProductTransformations trans;
00492
00493 if (trans.doTransform(sThis, sOther, rtn))
00494 {
00495 if (SymbolicTransformation::classVerbosity() > 0)
00496 {
00497 Out::println("Expr::operator*() transformed product\n["
00498 + toString() + "*"
00499 + other.toString() + "]\n to\n ["
00500 + rtn->toString() + "]");
00501 }
00502 return handle(rtn);
00503 }
00504
00505 return new ProductExpr(sThis, sOther);
00506 }
00507
00508 Expr Expr::operator-() const
00509 {
00510 TimeMonitor t(opTimer());
00511 TimeMonitor t1(unaryMinusTimer());
00512 Tabs tabs;
00513
00514 if (this->isComplex())
00515 {
00516 return new ComplexExpr(-real(), -imag());
00517 }
00518
00519
00520 if (this->size()==1)
00521 {
00522
00523 const SpectralExpr* se
00524 = dynamic_cast<const SpectralExpr*>((*this)[0].ptr().get());
00525 if (se != 0)
00526 {
00527 SpectralBasis basis = se->getSpectralBasis();
00528 Array<Expr> coeff(basis.nterms());
00529
00530 for(int i=0; i<basis.nterms(); i++)
00531 {
00532 coeff[i] = - se->getCoeff(i);
00533 }
00534 Expr rtn = new SpectralExpr( basis, coeff);
00535 return rtn;
00536 }
00537
00538
00539 const ConstantExpr* c
00540 = dynamic_cast<const ConstantExpr*>((*this)[0].ptr().get());
00541 const UnaryMinus* u
00542 = dynamic_cast<const UnaryMinus*>((*this)[0].ptr().get());
00543
00544 if (c != 0)
00545 {
00546 if (c->value()==0.0)
00547 {
00548 return new ZeroExpr();
00549 }
00550 else
00551 {
00552 return new ConstantExpr(-1.0 * c->value());
00553 }
00554 }
00555 else if (u != 0)
00556 {
00557 return u->arg();
00558 }
00559 else
00560 {
00561 RCP<ScalarExpr> sThis
00562 = rcp_dynamic_cast<ScalarExpr>((*this)[0].ptr());
00563 TEUCHOS_TEST_FOR_EXCEPTION(sThis.get()==NULL, std::logic_error,
00564 "Expr::operator-(): Operand " << (*this)[0].toString()
00565 << " is a non-scalar expression. All list structure "
00566 "should have been handled before this point");
00567 return new UnaryMinus(sThis);
00568 }
00569 }
00570
00571
00572 Array<Expr> rtn(this->size());
00573 for (int i=0; i<this->size(); i++)
00574 {
00575 rtn[i] = -((*this)[i]);
00576 }
00577 return new ListExpr(rtn);
00578 }
00579
00580
00581 Expr Expr::operator/(const Expr& other) const
00582 {
00583 TimeMonitor t(opTimer());
00584
00585
00586
00587 TEUCHOS_TEST_FOR_EXCEPTION(other.size() != 1,
00588 std::runtime_error,
00589 "Expr::operator/ detected division by a non-scalar "
00590 "expression " << toString());
00591
00592 TEUCHOS_TEST_FOR_EXCEPTION(other.isSpectral(), std::logic_error,
00593 "Division by a Spectral Expr is not yet defined");
00594
00595
00596 if (other.isComplex())
00597 {
00598 Expr magSq = other.real()*other.real() + other.imag()*other.imag();
00599 return (*this) * other.conj() / magSq;
00600 }
00601
00602
00603 if (isComplex() && !other.isComplex())
00604 {
00605 return new ComplexExpr(real()/other, imag()/other);
00606 }
00607
00608
00609 if (isSpectral() && !other.isSpectral() && !other.isComplex())
00610 {
00611 const SpectralExpr* se
00612 = dynamic_cast<const SpectralExpr*>((*this)[0].ptr().get());
00613
00614 SpectralBasis basis = se->getSpectralBasis();
00615 Array<Expr> coeff(basis.nterms());
00616
00617 for(int i=0; i<basis.nterms(); i++)
00618 {
00619 coeff[i] = se->getCoeff(i)/ other[0];
00620 }
00621 Expr rtn = new SpectralExpr( basis, coeff);
00622 return rtn;
00623 }
00624
00625
00626 if (this->size()==1)
00627 {
00628 return (*this)[0].divide(other[0]);
00629 }
00630
00631
00632 Array<Expr> rtn(this->size());
00633 for (int i=0; i<this->size(); i++)
00634 {
00635 rtn[i] = (*this)[i] / other;
00636 }
00637 return new ListExpr(rtn);
00638 }
00639
00640 const Expr& Expr::operator[](int i) const
00641 {
00642 TEUCHOS_TEST_FOR_EXCEPTION(ptr().get()==NULL, std::logic_error,
00643 "null this detected in Expr::operator[].");
00644
00645 TEUCHOS_TEST_FOR_EXCEPTION(i<0 || i>=(int) this->size(), std::runtime_error,
00646 "Expr::operator[]() index i=" << i << " out of range [0, "
00647 << this->size() << " in expr " << toString());
00648
00649 const ListExpr* le = dynamic_cast<const ListExpr*>(ptr().get());
00650
00651 if (le != 0)
00652 {
00653 const Expr& rtn = le->element(i);
00654 TEUCHOS_TEST_FOR_EXCEPTION(rtn.size() == 0, std::logic_error,
00655 "attempt to access an empty list; this should "
00656 "never happen because the case should have been "
00657 "dealt with earlier");
00658 if (rtn.size()==1)
00659 {
00660 TEUCHOS_TEST_FOR_EXCEPTION(rtn[0].ptr().get()==NULL, std::logic_error,
00661 "null return detected in Expr::operator[]. This="
00662 << toString() << ", i=" << i);
00663 return rtn[0];
00664 }
00665 TEUCHOS_TEST_FOR_EXCEPTION(rtn.ptr().get()==NULL, std::logic_error,
00666 "null return detected in Expr::operator[]. This="
00667 << toString() << ", i=" << i);
00668 return rtn;
00669 }
00670
00671 return *this;
00672 }
00673
00674 void Expr::append(const Expr& expr)
00675 {
00676 ListExpr* le = dynamic_cast<ListExpr*>(ptr().get());
00677
00678 if (le != 0)
00679 {
00680 le->append(expr);
00681 return;
00682 }
00683 else
00684 {
00685 if (ptr().get()==0)
00686 {
00687 Array<Expr> e(1);
00688 e[0] = expr;
00689 ptr() = rcp(new ListExpr(e));
00690 }
00691 else
00692 {
00693 Array<Expr> e(2);
00694 e[0] = *this;
00695 e[1] = expr;
00696 ptr() = rcp(new ListExpr(e));
00697 }
00698 }
00699 }
00700
00701 Expr Expr::flatten() const
00702 {
00703 const ListExpr* le = dynamic_cast<const ListExpr*>(ptr().get());
00704
00705 if (le != 0)
00706 {
00707 return le->flatten();
00708 }
00709 return *this;
00710 }
00711
00712 Expr Expr::flattenSpectral() const
00713 {
00714 Array<Expr> rtn(size());
00715 for (int i=0; i<size(); i++)
00716 {
00717 if ((*this)[i].size() == 1)
00718 {
00719 const SpectralExpr* se
00720 = dynamic_cast<const SpectralExpr*>((*this)[i][0].ptr().get());
00721 if (se != 0)
00722 {
00723 int nt = se->getSpectralBasis().nterms();
00724 Array<Expr> e(nt);
00725 for (int j=0; j<nt; j++)
00726 {
00727 e[j] = se->getCoeff(j);
00728 }
00729 rtn[i] = new ListExpr(e);
00730 }
00731 else
00732 {
00733 rtn[i] = (*this)[i];
00734 }
00735 }
00736 else
00737 {
00738 rtn[i] = (*this)[i].flattenSpectral();
00739 }
00740 }
00741 Expr r = new ListExpr(rtn);
00742 return r.flatten();
00743 }
00744
00745 int Expr::size() const
00746 {
00747 if (ptr().get()==0) return 0;
00748
00749 const ListExpr* le = dynamic_cast<const ListExpr*>(ptr().get());
00750
00751 if (le != 0)
00752 {
00753 return le->size();
00754 }
00755 return 1;
00756 }
00757
00758 int Expr::totalSize() const
00759 {
00760 if (ptr().get()==0) return 0;
00761
00762 const ListExpr* le = dynamic_cast<const ListExpr*>(ptr().get());
00763
00764 if (le != 0)
00765 {
00766 return le->totalSize();
00767 }
00768 return 1;
00769 }
00770
00771 Expr Expr::real() const
00772 {
00773 if (isComplex())
00774 {
00775 const ComplexExpr* cx = dynamic_cast<const ComplexExpr*>(ptr().get());
00776 return cx->real();
00777 }
00778 else if (size() > 1)
00779 {
00780 Array<Expr> rtn(size());
00781 for (int i=0; i<size(); i++)
00782 {
00783 rtn[i] = (*this)[i].real();
00784 }
00785 return new ListExpr(rtn);
00786 }
00787 else
00788 {
00789 return *this;
00790 }
00791 }
00792
00793 Expr Expr::imag() const
00794 {
00795 if (isComplex())
00796 {
00797 const ComplexExpr* cx = dynamic_cast<const ComplexExpr*>(ptr().get());
00798 return cx->imag();
00799 }
00800 else if (size() > 1)
00801 {
00802 Array<Expr> rtn(size());
00803 for (int i=0; i<size(); i++)
00804 {
00805 rtn[i] = (*this)[i].imag();
00806 }
00807 return new ListExpr(rtn);
00808 }
00809 else
00810 {
00811 return 0.0;
00812 }
00813
00814 }
00815
00816 Expr Expr::conj() const
00817 {
00818 if (size()==1)
00819 {
00820 if (isComplex())
00821 {
00822 return new ComplexExpr(real(), -imag());
00823 }
00824 else return real();
00825 }
00826 else
00827 {
00828 Array<Expr> rtn(size());
00829 for (int i=0; i<size(); i++)
00830 {
00831 rtn[i] = (*this)[i].conj();
00832 }
00833 return new ListExpr(rtn);
00834 }
00835 }
00836
00837 void Expr::setParameterValue(const double& value)
00838 {
00839 Parameter* pe = dynamic_cast<Parameter*>(ptr().get());
00840 TEUCHOS_TEST_FOR_EXCEPTION(pe==0, std::runtime_error,
00841 "Expr " << *this << " is not a Parameter expr, and "
00842 "so setParameterValue() should not be called");
00843 pe->setValue(value);
00844 }
00845
00846 double Expr::getParameterValue() const
00847 {
00848 const Parameter* pe = dynamic_cast<const Parameter*>(ptr().get());
00849 TEUCHOS_TEST_FOR_EXCEPTION(pe==0, std::runtime_error,
00850 "Expr " << *this << " is not a Parameter expr, and "
00851 "so getParameterValue() should not be called");
00852 return pe->value();
00853 }
00854
00855
00856 bool Expr::isIndependentOf(const Expr& u) const
00857 {
00858 TEUCHOS_TEST_FOR_EXCEPTION(ptr().get()==0, std::runtime_error,
00859 "function called on null expression");
00860
00861 return scalarExpr()->isIndependentOf(u);
00862 }
00863
00864
00865 bool Expr::isLinearForm(const Expr& u) const
00866 {
00867 TEUCHOS_TEST_FOR_EXCEPTION(ptr().get()==0, std::runtime_error,
00868 "function called on null expression");
00869
00870 return scalarExpr()->isLinearForm(u);
00871 }
00872
00873
00874
00875 bool Expr::isQuadraticForm(const Expr& u) const
00876 {
00877 TEUCHOS_TEST_FOR_EXCEPTION(ptr().get()==0, std::runtime_error,
00878 "function called on null expression");
00879
00880 return scalarExpr()->isQuadraticForm(u);
00881 }
00882
00883
00884
00885 bool Expr::isLinearInTests() const
00886 {
00887 TEUCHOS_TEST_FOR_EXCEPTION(ptr().get()==0, std::runtime_error,
00888 "function called on null expression");
00889
00890 return scalarExpr()->isLinearInTests();
00891 }
00892
00893
00894 bool Expr::hasTestFunctions() const
00895 {
00896 TEUCHOS_TEST_FOR_EXCEPTION(ptr().get()==0, std::runtime_error,
00897 "function called on null expression");
00898
00899 return scalarExpr()->hasTestFunctions();
00900 }
00901
00902
00903 bool Expr::everyTermHasTestFunctions() const
00904 {
00905 TEUCHOS_TEST_FOR_EXCEPTION(ptr().get()==0, std::runtime_error,
00906 "function called on null expression");
00907
00908 return scalarExpr()->everyTermHasTestFunctions();
00909 }
00910
00911 bool Expr::isTestElement() const
00912 {
00913 const TestFuncElement* uPtr
00914 = dynamic_cast<const TestFuncElement*>(ptr().get());
00915 return uPtr != 0;
00916 }
00917
00918
00919 bool Expr::isUnknownElement() const
00920 {
00921 const UnknownFuncElement* uPtr
00922 = dynamic_cast<const UnknownFuncElement*>(ptr().get());
00923 return uPtr != 0;
00924 }
00925
00926
00927 void Expr::getUnknowns(Set<int>& unkID,
00928 Array<Expr>& unks) const
00929 {
00930 TEUCHOS_TEST_FOR_EXCEPTION(ptr().get()==0, std::runtime_error,
00931 "function called on null expression");
00932
00933 const UnknownFuncElement* u
00934 = dynamic_cast<const UnknownFuncElement*>(ptr().get());
00935
00936 if (u != 0)
00937 {
00938 if (!unkID.contains(u->fid().dofID()))
00939 {
00940 unks.append(*this);
00941 unkID.put(u->fid().dofID());
00942 }
00943 }
00944 else
00945 {
00946 scalarExpr()->getUnknowns(unkID, unks);
00947 }
00948 }
00949
00950
00951 void Expr::getTests(Set<int>& varID, Array<Expr>& vars) const
00952 {
00953 TEUCHOS_TEST_FOR_EXCEPTION(ptr().get()==0, std::runtime_error,
00954 "function called on null expression");
00955
00956 const TestFuncElement* u
00957 = dynamic_cast<const TestFuncElement*>(ptr().get());
00958
00959 if (u != 0)
00960 {
00961 if (!varID.contains(u->fid().dofID()))
00962 {
00963 vars.append(*this);
00964 varID.put(u->fid().dofID());
00965 }
00966 }
00967 else
00968 {
00969 scalarExpr()->getTests(varID, vars);
00970 }
00971 }
00972
00973
00974 const ScalarExpr* Expr::scalarExpr() const
00975 {
00976 const ScalarExpr* se = dynamic_cast<const ScalarExpr*>(ptr().get());
00977 TEUCHOS_TEST_FOR_EXCEPTION(se==0, std::logic_error, "ScalarExpr expected here");
00978 return se;
00979 }
00980
00981 const FuncElementBase* Expr::funcElement() const
00982 {
00983 const FuncElementBase* fe
00984 = dynamic_cast<const FuncElementBase*>(ptr().get());
00985 return fe;
00986 }
00987
00988
00989
00990
00991 namespace Sundance
00992 {
00993 using namespace Sundance;
00994
00995 Expr Complex(const Expr& re, const Expr& im)
00996 {
00997 TEUCHOS_TEST_FOR_EXCEPTION(re.size() != im.size(), std::runtime_error,
00998 "arguments mismatched in Complex(). Real part="
00999 << re << ", imaginary part=" << im);
01000
01001 TEUCHOS_TEST_FOR_EXCEPTION(re.isComplex() || im.isComplex(), std::runtime_error,
01002 "recursively defined complex number. Real part="
01003 << re << ", imaginary part=" << im);
01004
01005 if (re.totalSize() > 1)
01006 {
01007 Array<Expr> rtn(re.size());
01008 for (int i=0; i<re.size(); i++)
01009 {
01010 rtn[i] = Complex(re[i], im[i]);
01011 }
01012 return new ListExpr(rtn);
01013 }
01014
01015 const ZeroExpr* zr = dynamic_cast<const ZeroExpr*>(re[0].ptr().get());
01016 const ZeroExpr* zi = dynamic_cast<const ZeroExpr*>(im[0].ptr().get());
01017
01018 if (zr == 0)
01019 {
01020 if (zi==0)
01021 {
01022 return new ComplexExpr(re, im);
01023 }
01024 else
01025 {
01026 return re;
01027 }
01028 }
01029 else
01030 {
01031 if (zi != 0)
01032 {
01033 return Expr(0.0);
01034 }
01035 else
01036 {
01037 return new ComplexExpr(0.0, im);
01038 }
01039 }
01040 return new ComplexExpr(re, im);
01041 }
01042
01043 Expr List(const Expr& a)
01044 {
01045 return new ListExpr(tuple(a));
01046 }
01047
01048 Expr List(const Expr& a, const Expr& b)
01049 {
01050 return new ListExpr(tuple(a,b));
01051 }
01052
01053 Expr List(const Expr& a, const Expr& b, const Expr& c)
01054 {
01055 return new ListExpr(tuple(a,b,c));
01056 }
01057
01058 Expr List(const Expr& a, const Expr& b, const Expr& c,
01059 const Expr& d)
01060 {
01061 return new ListExpr(tuple(a,b,c,d));
01062 }
01063
01064 Expr List(const Expr& a, const Expr& b, const Expr& c,
01065 const Expr& d, const Expr& e)
01066 {
01067 return new ListExpr(tuple(a,b,c,d,e));
01068 }
01069
01070 Expr List(const Expr& a, const Expr& b, const Expr& c,
01071 const Expr& d, const Expr& e, const Expr& f)
01072 {
01073 return new ListExpr(tuple(a,b,c,d,e,f));
01074 }
01075
01076 Expr List(const Expr& a, const Expr& b, const Expr& c,
01077 const Expr& d, const Expr& e, const Expr& f,
01078 const Expr& g)
01079 {
01080 return new ListExpr(tuple(a,b,c,d,e,f,g));
01081 }
01082
01083 Expr List(const Expr& a, const Expr& b, const Expr& c,
01084 const Expr& d, const Expr& e, const Expr& f,
01085 const Expr& g, const Expr& h)
01086 {
01087 return new ListExpr(tuple(a,b,c,d,e,f,g,h));
01088 }
01089
01090
01091 }
01092
