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 "SundanceDiffOp.hpp"
00032 #include "PlayaExceptions.hpp"
00033
00034 #include "SundanceTestFuncElement.hpp"
00035
00036 #include "SundanceCoordExpr.hpp"
00037 #include "SundanceZeroExpr.hpp"
00038 #include "PlayaTabs.hpp"
00039 #include "SundanceOut.hpp"
00040
00041 using namespace Sundance;
00042 using namespace Sundance;
00043
00044 using namespace Sundance;
00045 using namespace Teuchos;
00046
00047
00048 DiffOp::DiffOp(const MultiIndex& op,
00049 const RCP<ScalarExpr>& arg)
00050 : UnaryExpr(arg), mi_(op), myCoordDeriv_(coordDeriv(op.firstOrderDirection())), requiredFunctions_(),
00051 ignoreFuncTerms_(false)
00052 {}
00053
00054 void DiffOp::registerSpatialDerivs(const EvalContext& context,
00055 const Set<MultiIndex>& miSet) const
00056 {
00057 EvaluatableExpr::registerSpatialDerivs(context, miSet);
00058 Set<MultiIndex> s;
00059 for (Set<MultiIndex>::const_iterator i=miSet.begin(); i!=miSet.end(); i++)
00060 {
00061 const MultiIndex& m = *i;
00062 s.put(m+mi_);
00063 }
00064 evaluatableArg()->registerSpatialDerivs(context, s);
00065 }
00066
00067
00068 std::ostream& DiffOp::toText(std::ostream& os, bool ) const
00069 {
00070 std::string miStr = CoordExpr::coordName(mi_.firstOrderDirection(), "");
00071 os << "D[" << arg().toString() << ", " << miStr << "]";
00072 return os;
00073 }
00074
00075 XMLObject DiffOp::toXML() const
00076 {
00077 XMLObject rtn("DiffOp");
00078 rtn.addAttribute("m", mi_.toString());
00079 rtn.addChild(arg().toXML());
00080
00081 return rtn;
00082 }
00083
00084
00085 Evaluator* DiffOp::createEvaluator(const EvaluatableExpr* expr,
00086 const EvalContext& context) const
00087 {
00088 return new DiffOpEvaluator(dynamic_cast<const DiffOp*>(expr), context);
00089 }
00090
00091
00092
00093 void DiffOp::requestMultiIndexAtEvalPoint(const MultiIndex& mi,
00094 const MultipleDeriv& u,
00095 const EvalContext& context) const
00096 {
00097 int verb = context.setupVerbosity();
00098 Tabs tab0(0);
00099 SUNDANCE_MSG3(verb, tab0 << "DiffOp::requestMultiIndexAtEvalPoint() for=" << toString());
00100 TEUCHOS_TEST_FOR_EXCEPT(u.size() != 1);
00101 const Deriv& d = *(u.begin());
00102
00103 if (d.isFunctionalDeriv())
00104 {
00105 const SpatialDerivSpecifier& sds = d.opOnFunc();
00106
00107 TEUCHOS_TEST_FOR_EXCEPTION(sds.isDivergence(), std::logic_error,
00108 "divergence operators not possible within DiffOp");
00109 TEUCHOS_TEST_FOR_EXCEPTION(sds.isNormal(), std::logic_error,
00110 "normal deriv operators not possible within DiffOp");
00111
00112 const MultiIndex& newMi = sds.mi();
00113
00114 const SymbolicFuncElement* sfe = d.symbFuncElem();
00115 TEUCHOS_TEST_FOR_EXCEPT(sfe == 0);
00116 const EvaluatableExpr* evalPt = sfe->evalPt();
00117 const ZeroExpr* z = dynamic_cast<const ZeroExpr*>(evalPt);
00118 if (z != 0) return;
00119 const DiscreteFuncElement* df
00120 = dynamic_cast<const DiscreteFuncElement*>(evalPt);
00121 df->addMultiIndex(newMi);
00122 df->findW(1, context);
00123 df->findV(1, context);
00124 df->findC(1, context);
00125 }
00126 }
00127
00128
00129 RCP<Array<Set<MultipleDeriv> > >
00130 DiffOp::internalDetermineR(const EvalContext& context,
00131 const Array<Set<MultipleDeriv> >& RInput) const
00132 {
00133 Tabs tab0(0);
00134 int verb = context.setupVerbosity();
00135
00136 SUNDANCE_MSG2(verb, tab0 << "DiffOp::internalDetermineR for=" << toString());
00137 SUNDANCE_MSG2(verb, tab0 << "RInput = " << RInput );
00138
00139 RCP<Array<Set<MultipleDeriv> > > rtn
00140 = rcp(new Array<Set<MultipleDeriv> >(RInput.size()));
00141
00142 {
00143 Tabs tab1;
00144 for (int i=0; i<RInput.size(); i++)
00145 {
00146 Tabs tab2;
00147 const Set<MultipleDeriv>& Wi = findW(i, context);
00148 SUNDANCE_MSG5(verb, tab2 << "W[" << i << "] = " << Wi );
00149 (*rtn)[i] = RInput[i].intersection(Wi);
00150 }
00151
00152 const Set<MultipleDeriv>& W1 = evaluatableArg()->findW(1, context);
00153 SUNDANCE_MSG3(verb, tab1 << "arg W1 = " << W1);
00154 Set<MultipleDeriv> ZxXx = applyZx(W1, mi_).setUnion(Xx(mi_));
00155 SUNDANCE_MSG3(verb, tab1 << "Z union X = " << ZxXx);
00156 Array<Set<MultipleDeriv> > RArg(RInput.size()+1);
00157 RArg[0].put(MultipleDeriv());
00158 RArg[1].put(MultipleDeriv(coordDeriv(mi_.firstOrderDirection())));
00159
00160
00161
00162 for (int order=0; order<RInput.size(); order++)
00163 {
00164 Tabs tab2;
00165 const Set<MultipleDeriv>& WArgPlus = evaluatableArg()->findW(order+1, context);
00166 const Set<MultipleDeriv>& WArg = evaluatableArg()->findW(order, context);
00167 SUNDANCE_MSG3(verb, tab2 << "order = " << order);
00168 SUNDANCE_MSG3(verb, tab2 << "RInput = " << RInput[order]);
00169 SUNDANCE_MSG3(verb, tab2 << "WArg = " << WArg);
00170 SUNDANCE_MSG3(verb, tab2 << "WArgPlus = " << WArgPlus);
00171 SUNDANCE_MSG3(verb, tab2 << "ZxXx times RInput = "
00172 << setProduct(ZxXx, RInput[order]));
00173 SUNDANCE_MSG3(verb, tab2 << "Tx(RInput, " << -mi_ << ") = "
00174 << applyTx(RInput[order], -mi_) );
00175 RArg[order+1].merge(setProduct(ZxXx, RInput[order]).intersection(WArgPlus));
00176 RArg[order].merge(applyTx(RInput[order], -mi_).intersection(WArg));
00177 }
00178 SUNDANCE_MSG3(verb, tab1 << "RArg = " << RArg);
00179
00180 SUNDANCE_MSG2(verb, tab1 << "calling determineR() for arg "
00181 << evaluatableArg()->toString());
00182 evaluatableArg()->determineR(context, RArg);
00183
00184
00185
00186
00187 SUNDANCE_MSG3(verb, tab1 << "calling findR(1) to determine required spatial derivatives ");
00188 const Set<MultipleDeriv>& RArg1 = evaluatableArg()->findR(1, context);
00189 SUNDANCE_MSG3(verb, tab1 << "RArg1 = " << RArg1);
00190 for (Set<MultipleDeriv>::const_iterator i=RArg1.begin(); i!=RArg1.end(); i++)
00191 {
00192 requestMultiIndexAtEvalPoint(mi(), *i, context);
00193 }
00194 }
00195 printR(verb, rtn);
00196 SUNDANCE_MSG2(verb, tab0 << "done with DiffOp::internalDetermineR for "
00197 << toString());
00198
00199 return rtn;
00200 }
00201
00202
00203 Set<MultipleDeriv> DiffOp::internalFindW(int order, const EvalContext& context) const
00204 {
00205 Tabs tabs(0);
00206 int verb = context.setupVerbosity();
00207
00208 SUNDANCE_MSG2(verb, tabs << "DiffOp::internalFindW(order="
00209 << order << ") for " << toString());
00210
00211 Set<MultipleDeriv> rtn ;
00212 if (order <= 2)
00213 {
00214 Tabs tab1;
00215 const Set<MultipleDeriv>& W1 = evaluatableArg()->findW(1, context);
00216 const Set<MultipleDeriv>& WArg = evaluatableArg()->findW(order, context);
00217 const Set<MultipleDeriv>& WArgPlus = evaluatableArg()->findW(order+1, context);
00218
00219 SUNDANCE_MSG5(verb, tab1 << "W1 = " << W1);
00220 SUNDANCE_MSG5(verb, tab1 << "WArg = " << WArg);
00221 SUNDANCE_MSG5(verb, tab1 << "WArgPlus = " << WArgPlus);
00222 Set<MultipleDeriv> Tx = applyTx(WArg, mi_);
00223 Set<MultipleDeriv> ZXx = applyZx(W1, mi_).setUnion(Xx(mi_));
00224 SUNDANCE_MSG5(verb, tab1 << "Tx(Warg) = " << Tx);
00225 SUNDANCE_MSG5(verb, tab1 << "ZXx(W1) = " << ZXx);
00226 Set<MultipleDeriv> WargPlusOslashZXx = setDivision(WArgPlus, ZXx);
00227 SUNDANCE_MSG5(verb, tab1 << "WArgPlus / ZXx = "
00228 << WargPlusOslashZXx);
00229 rtn = WargPlusOslashZXx.setUnion(Tx);
00230 }
00231 SUNDANCE_MSG3(verb, tabs << "W[" << order << "]=" << rtn);
00232 SUNDANCE_MSG3(verb, tabs << "done with DiffOp::internalFindW(" << order << ") for "
00233 << toString());
00234
00235 return rtn;
00236 }
00237
00238
00239 Set<MultipleDeriv> DiffOp::internalFindV(int order, const EvalContext& context) const
00240 {
00241 Tabs tabs(0);
00242 int verb = context.setupVerbosity();
00243 SUNDANCE_MSG2(verb, tabs << "DiffOp::internalFindV(" << order << ") for "
00244 << toString());
00245
00246 Set<MultipleDeriv> rtn;
00247
00248 if (order <= 2)
00249 {
00250 Tabs tab1;
00251 const Set<MultipleDeriv>& W1 = evaluatableArg()->findW(1, context);
00252 const Set<MultipleDeriv>& VArg = evaluatableArg()->findV(order, context);
00253 const Set<MultipleDeriv>& VArgPlus
00254 = evaluatableArg()->findV(order+1, context);
00255 const Set<MultipleDeriv>& WArg = evaluatableArg()->findW(order, context);
00256 const Set<MultipleDeriv>& WArgPlus
00257 = evaluatableArg()->findW(order+1, context);
00258
00259 SUNDANCE_MSG5(verb, tab1 << "W1=" << W1);
00260 SUNDANCE_MSG5(verb, tab1 << "VArg=" << VArg);
00261 SUNDANCE_MSG5(verb, tab1 << "VArgPlus=" << VArgPlus);
00262 SUNDANCE_MSG5(verb, tab1 << "WArg=" << WArg);
00263 SUNDANCE_MSG5(verb, tab1 << "WArgPlus=" << WArgPlus);
00264
00265 Set<MultipleDeriv> Tx = applyTx(VArg, mi_);
00266 Set<MultipleDeriv> Zx = applyZx(W1, mi_);
00267 SUNDANCE_MSG5(verb, tab1 << "Tx(Varg) = " << Tx);
00268 SUNDANCE_MSG5(verb, tab1 << "Zx(W1) = " << Zx);
00269 Set<MultipleDeriv> WargPlusOslashZx = setDivision(WArgPlus, Zx);
00270 Set<MultipleDeriv> VargPlusOslashXx = setDivision(VArgPlus, Xx(mi_));
00271 SUNDANCE_MSG5(verb, tab1 << "WArgPlus / Z_alpha = "
00272 << WargPlusOslashZx);
00273 SUNDANCE_MSG5(verb, tab1 << "VArgPlus / X_alpha = "
00274 << VargPlusOslashXx);
00275 rtn = WargPlusOslashZx.setUnion(VargPlusOslashXx).setUnion(Tx);
00276
00277 SUNDANCE_MSG5(verb, tab1 << "WArgPlus/Z union VArgPlus/X union T =" << rtn);
00278 rtn = rtn.intersection(findR(order, context));
00279 }
00280
00281 SUNDANCE_MSG2(verb, tabs << "V[" << order << "]=" << rtn);
00282 SUNDANCE_MSG2(verb, tabs << "done with DiffOp::internalFindV(" << order << ") for "
00283 << toString());
00284
00285 return rtn;
00286 }
00287
00288
00289 Set<MultipleDeriv> DiffOp::internalFindC(int order, const EvalContext& context) const
00290 {
00291 Tabs tabs(0);
00292 int verb = context.setupVerbosity();
00293 SUNDANCE_MSG2(verb, tabs << "DiffOp::internalFindC() for "
00294 << toString());
00295 Set<MultipleDeriv> rtn ;
00296
00297 {
00298 Tabs tab1;
00299 SUNDANCE_MSG5(verb, tab1 << "finding R");
00300 const Set<MultipleDeriv>& R = findR(order, context);
00301 SUNDANCE_MSG5(verb, tab1 << "finding V");
00302 const Set<MultipleDeriv>& V = findV(order, context);
00303
00304 evaluatableArg()->findC(order, context);
00305
00306 SUNDANCE_MSG5(verb, tab1 << "R=" << R);
00307 SUNDANCE_MSG5(verb, tab1 << "V=" << V);
00308 rtn = R.setDifference(V);
00309 SUNDANCE_MSG3(verb, tabs << "C[" << order << "]=" << rtn);
00310 }
00311
00312 SUNDANCE_MSG2(verb, tabs << "C[" << order << "]=R\\V = " << rtn);
00313 SUNDANCE_MSG2(verb, tabs << "done with DiffOp::internalFindC for "
00314 << toString());
00315 return rtn;
00316 }
00317
00318
00319
00320
00321 bool DiffOp::lessThan(const ScalarExpr* other) const
00322 {
00323 const DiffOp* d = dynamic_cast<const DiffOp*>(other);
00324 TEUCHOS_TEST_FOR_EXCEPTION(d==0, std::logic_error, "cast should never fail at this point");
00325
00326 if (myCoordDeriv_ < d->myCoordDeriv_) return true;
00327 if (d->myCoordDeriv_ < myCoordDeriv_) return false;
00328
00329 return ExprWithChildren::lessThan(other);
00330 }
00331