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 "SundanceMultipleDeriv.hpp"
00032 #include "SundanceSymbolicFuncElement.hpp"
00033
00034 using namespace Sundance;
00035 using namespace Sundance;
00036 using namespace Teuchos;
00037
00038 MultipleDeriv::MultipleDeriv()
00039 : MultiSet<Deriv>()
00040 {}
00041
00042 MultipleDeriv::MultipleDeriv(const Deriv& d)
00043 : MultiSet<Deriv>()
00044 {
00045 put(d);
00046 }
00047 MultipleDeriv::MultipleDeriv(const Deriv& d1, const Deriv& d2)
00048 : MultiSet<Deriv>()
00049 {
00050 put(d1);
00051 put(d2);
00052 }
00053
00054 int MultipleDeriv::spatialOrder() const
00055 {
00056 int rtn = 0;
00057 for (MultipleDeriv::const_iterator i=this->begin(); i!=this->end(); i++)
00058 {
00059 if (i->isCoordDeriv())
00060 {
00061 rtn += 1;
00062 }
00063 }
00064 return rtn;
00065 }
00066
00067 MultiIndex MultipleDeriv::spatialDeriv() const
00068 {
00069 MultiIndex rtn;
00070 for (MultipleDeriv::const_iterator i=this->begin(); i!=this->end(); i++)
00071 {
00072 if (i->isCoordDeriv())
00073 {
00074 int d = i->coordDerivDir();
00075 rtn[d] += 1;
00076 }
00077 }
00078 return rtn;
00079 }
00080
00081 MultiSet<FunctionIdentifier> MultipleDeriv::funcIDs() const
00082 {
00083 MultiSet<FunctionIdentifier> rtn;
00084 for (MultipleDeriv::const_iterator i=this->begin(); i!=this->end(); i++)
00085 {
00086 if (i->isFunctionalDeriv())
00087 {
00088 rtn.put(i->fid());
00089 }
00090 TEUCHOS_TEST_FOR_EXCEPTION(!i->isFunctionalDeriv(), std::runtime_error,
00091 "MultipleDeriv::funcIDs() found spatial deriv");
00092 }
00093 return rtn;
00094 }
00095
00096 MultiSet<int> MultipleDeriv::dofIDs() const
00097 {
00098 MultiSet<int> rtn;
00099 for (MultipleDeriv::const_iterator i=this->begin(); i!=this->end(); i++)
00100 {
00101 if (i->isFunctionalDeriv())
00102 {
00103 int f = i->dofID();
00104 rtn.put(f);
00105 }
00106 TEUCHOS_TEST_FOR_EXCEPTION(!i->isFunctionalDeriv(), std::runtime_error,
00107 "MultipleDeriv::sharedFuncIDs() found spatial deriv");
00108 }
00109 return rtn;
00110 }
00111
00112 MultipleDeriv MultipleDeriv::product(const MultipleDeriv& other) const
00113 {
00114 MultipleDeriv rtn;
00115
00116 for (MultipleDeriv::const_iterator i=this->begin(); i!=this->end(); i++)
00117 {
00118 rtn.put(*i);
00119 }
00120 for (MultipleDeriv::const_iterator i=other.begin(); i!=other.end(); i++)
00121 {
00122 rtn.put(*i);
00123 }
00124 return rtn;
00125 }
00126
00127 MultipleDeriv MultipleDeriv::factorOutDeriv(const Deriv& x) const
00128 {
00129 MultipleDeriv rtn = *this;
00130
00131 MultipleDeriv::iterator i = rtn.find(x);
00132
00133
00134 if (i != rtn.end()) rtn.erase(i);
00135
00136 if (rtn.size() == this->size()) return MultipleDeriv();
00137
00138 return rtn;
00139 }
00140
00141
00142 bool MultipleDeriv::containsDeriv(const MultipleDeriv& x) const
00143 {
00144 for (MultipleDeriv::const_iterator i=x.begin(); i!=x.end(); i++)
00145 {
00146 if ( count(*i) <= x.count(*i) ) return false;
00147 }
00148 return true;
00149 }
00150
00151 MultipleDeriv MultipleDeriv::factorOutDeriv(const MultipleDeriv& x) const
00152 {
00153 MultipleDeriv rtn = *this;
00154
00155 for (MultipleDeriv::const_iterator i=x.begin(); i!=x.end(); i++)
00156 {
00157 MultipleDeriv::iterator j = rtn.find(*i);
00158
00159
00160 if (j != rtn.end()) rtn.erase(j);
00161 }
00162
00163 if (rtn.size() == this->size()) return MultipleDeriv();
00164 return rtn;
00165 }
00166
00167 bool MultipleDeriv
00168 ::isInRequiredSet(const Set<MultiSet<int> >& funcCombinations,
00169 const Set<MultiIndex>& multiIndices) const
00170 {
00171 if (spatialOrder() == 0)
00172 {
00173 return funcCombinations.contains(dofIDs());
00174 }
00175 else
00176 {
00177 return multiIndices.contains(spatialDeriv());
00178 }
00179 }
00180
00181
00182 void MultipleDeriv
00183 ::productRulePermutations(ProductRulePerms& perms) const
00184 {
00185 int N = order();
00186
00187 if (N==0)
00188 {
00189 MultipleDeriv md0;
00190 DerivPair p(md0, md0);
00191 perms.put(p, 1);
00192 return;
00193 }
00194
00195 int p2 = pow2(N);
00196
00197 for (int i=0; i<p2; i++)
00198 {
00199 MultipleDeriv left;
00200 MultipleDeriv right;
00201 Array<int> bits = bitsOfAnInteger(i, N);
00202 int j=0;
00203 MultipleDeriv::const_iterator iter;
00204 for (iter=this->begin(); iter != this->end(); iter++, j++)
00205 {
00206 if (bits[j]==true)
00207 {
00208 left.put(*iter);
00209 }
00210 else
00211 {
00212 right.put(*iter);
00213 }
00214 }
00215 DerivPair p(left, right);
00216 if (!perms.containsKey(p))
00217 {
00218 perms.put(p, 1);
00219 }
00220 else
00221 {
00222 int count = perms.get(p);
00223 perms.put(p, count+1);
00224 }
00225 }
00226 }
00227
00228 Array<int> MultipleDeriv::bitsOfAnInteger(int x, int n)
00229 {
00230 TEUCHOS_TEST_FOR_EXCEPTION(x >= pow2(n), std::logic_error,
00231 "Invalid input to MultipleDeriv::bitsOfX");
00232
00233 Array<int> rtn(n);
00234
00235 int r = x;
00236 for (int b=n-1; b>=0; b--)
00237 {
00238 rtn[b] = r/pow2(b);
00239 r = r - rtn[b]*pow2(b);
00240 }
00241 return rtn;
00242 }
00243
00244 int MultipleDeriv::pow2(int n)
00245 {
00246 static Array<int> p2(1,1);
00247
00248 if (n >= p2.size())
00249 {
00250 int oldN = p2.size();
00251 for (int i=oldN; i<=n; i++) p2.push_back(p2[i-1]*2);
00252 }
00253
00254 return p2[n];
00255 }
00256
00257
00258
00259 namespace Sundance
00260 {
00261 Set<MultipleDeriv> applyTx(const Set<MultipleDeriv>& s,
00262 const MultiIndex& x)
00263 {
00264 Set<MultipleDeriv> rtn;
00265
00266 for (Set<MultipleDeriv>::const_iterator i=s.begin(); i!=s.end(); i++)
00267 {
00268 const MultipleDeriv& md = *i;
00269 for (MultipleDeriv::const_iterator j=md.begin(); j!=md.end(); j++)
00270 {
00271 const Deriv& d = *j;
00272 if (d.isFunctionalDeriv())
00273 {
00274 const MultiIndex& mi = d.opOnFunc().mi();
00275 MultiIndex miNew = mi+x;
00276 if (miNew.isValid())
00277 {
00278 Deriv dNew = d.derivWrtMultiIndex(miNew);
00279 MultipleDeriv mdNew = md;
00280 mdNew.erase(mdNew.find(d));
00281 mdNew.put(dNew);
00282 rtn.put(mdNew);
00283 }
00284 }
00285 }
00286 }
00287 return rtn;
00288 }
00289
00290 Set<MultipleDeriv> Xx(const MultiIndex& x)
00291 {
00292 Set<MultipleDeriv> rtn;
00293
00294 TEUCHOS_TEST_FOR_EXCEPTION(x.order() < 0 || x.order() > 1, std::logic_error,
00295 "invalid multiindex " << x << " in this context");
00296
00297 MultipleDeriv xmd = makeMultiDeriv(coordDeriv(x.firstOrderDirection()));
00298 rtn.put(xmd);
00299 return rtn;
00300 }
00301
00302 Set<MultipleDeriv> applyZx(const Set<MultipleDeriv>& W,
00303 const MultiIndex& x)
00304 {
00305 Set<MultipleDeriv> rtn;
00306
00307 TEUCHOS_TEST_FOR_EXCEPTION(x.order() < 0 || x.order() > 1, std::logic_error,
00308 "invalid multiindex " << x << " in this context");
00309
00310 for (Set<MultipleDeriv>::const_iterator i=W.begin(); i!=W.end(); i++)
00311 {
00312 const MultipleDeriv& md = *i;
00313 TEUCHOS_TEST_FOR_EXCEPTION(md.order() != 1, std::logic_error,
00314 "Only first-order multiple functional derivatives "
00315 "should appear in this function. The derivative "
00316 << md << " is not first-order.");
00317
00318 const Deriv& d = *(md.begin());
00319
00320 if (d.isFunctionalDeriv())
00321 {
00322
00323 TEUCHOS_TEST_FOR_EXCEPTION(!d.canBeDifferentiated(),
00324 std::logic_error, "function signature " << d << " cannot be "
00325 "differentiated further spatially");
00326
00327
00328 const SymbolicFuncElement* sfe = d.symbFuncElem();
00329 TEUCHOS_TEST_FOR_EXCEPTION(sfe==0, std::logic_error,
00330 "can't cast function in "
00331 << d << " to a SymbolicFuncElement");
00332 if (sfe && !sfe->evalPtIsZero()) rtn.put(md);
00333 }
00334 }
00335 return rtn;
00336 }
00337
00338
00339
00340 int factorial(const MultipleDeriv& ms)
00341 {
00342 Sundance::Map<Deriv, int> counts;
00343
00344 for (MultipleDeriv::const_iterator i=ms.begin(); i!=ms.end(); i++)
00345 {
00346 if (counts.containsKey(*i)) counts[*i]++;
00347 else counts.put(*i, 1);
00348 }
00349
00350 int rtn = 1;
00351 for (Sundance::Map<Deriv, int>::const_iterator
00352 i=counts.begin(); i!=counts.end(); i++)
00353 {
00354 int f = 1;
00355 for (int j=1; j<=i->second; j++) f *= j;
00356 rtn *= f;
00357 }
00358 return rtn;
00359 }
00360
00361 MultipleDeriv makeMultiDeriv(const Deriv& d)
00362 {
00363 MultipleDeriv rtn;
00364 rtn.put(d);
00365 return rtn;
00366 }
00367
00368 bool hasParameter(const MultipleDeriv& d)
00369 {
00370 for (MultipleDeriv::const_iterator i=d.begin(); i!=d.end(); i++)
00371 {
00372 if (i->isParameter()) return true;
00373 }
00374 return false;
00375 }
00376
00377
00378
00379 }