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 "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 /* paren */) 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     /* inform the evaluation points of all functions appearing in the argument
00185      * that we'll need their spatial derivatives in direction mi(). */
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   /* all done */  
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     /** Call findC() to ensure that the argument has C tabulated */
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
SundanceDiffOp.cpp
