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 #ifndef SUNDANCE_MULTIPLEDERIV_H 00032 #define SUNDANCE_MULTIPLEDERIV_H 00033 00034 #include "SundanceDefs.hpp" 00035 #include "SundanceDeriv.hpp" 00036 #include "SundanceMultiSet.hpp" 00037 #include "SundanceOrderedTuple.hpp" 00038 #include "SundanceMap.hpp" 00039 #include "Teuchos_Array.hpp" 00040 00041 00042 00043 namespace Sundance 00044 { 00045 using namespace Sundance; 00046 class MultipleDeriv; 00047 /** */ 00048 typedef OrderedPair<MultipleDeriv, MultipleDeriv> DerivPair; 00049 00050 /** */ 00051 typedef Sundance::Map<DerivPair, int> ProductRulePerms; 00052 /** Class MultipleDeriv is a multiple functional derivative operator 00053 * represented as a multiset of first-order derivatives. 00054 * The derivatives are each Deriv objects, so the multiple 00055 * derivative can be an arbitrary mixture of spatial and 00056 * functional derivatives. 00057 * 00058 * <h3> Product rule application </h3> 00059 * 00060 * Class MultipleDeriv contains a utility method for computing the 00061 * derivatives that appear in the general product rule, as 00062 * used in automatic differentiation of products and spatial 00063 * derivatives. 00064 * The arbitrary order, multivariable product rule gives 00065 * the derivative of a product \f$L\cdot R\f$ 00066 * with respect to a multiset of variables 00067 * \f$\{u_1, u_2, \dots, u_N\}\f$ as the sum 00068 * \f[ D_{u_1} D_{u_2} \cdots D_{u_N} (L \cdot R) 00069 * =\sum_{i=1}^{2^N} \left[\prod_{j=1}^N D_{u_j}^{b_{i,j}}\right] L 00070 * \times \left[\prod_{j=1}^N D_{u_j}^{1-b_{i,j}}\right] R\f] 00071 * where \f$b_{i,j}\f$ is the \f$j\f$-th bit in the binary 00072 * representation of \f$i\f$. Method productRulePermutations() 00073 * enumerates all the permutations that apply to the left and right 00074 * operands in this sum, and returns the 00075 * arrays of left operators and right operators 00076 * through non-const reference arguments. 00077 * The left and right arguments return 00078 * 00079 * \f[ {\rm left = Array}_{i=1}^{2^N} 00080 * \{\prod_{j=1}^N D_{u_j}^{b_{i,j}}\}\f] 00081 * 00082 * \f[ {\rm right = Array}_{i=1}^{2^N} 00083 * \{\prod_{j=1}^N D_{u_j}^{1-b_{i,j}}\}\f] 00084 * 00085 * Each element in these arrays is a multiple derivative, and 00086 * is naturally represented with a MultipleDeriv object. 00087 */ 00088 class MultipleDeriv : public MultiSet<Deriv> 00089 { 00090 public: 00091 00092 /** Construct an empty multiple derivative */ 00093 MultipleDeriv(); 00094 00095 /** Construct a first-order multiple deriv */ 00096 MultipleDeriv(const Deriv& d); 00097 00098 /** Construct a second-order multiple deriv */ 00099 MultipleDeriv(const Deriv& d1, const Deriv& d2); 00100 00101 /** Return the order of this derivative. Since a 00102 * multiple derivative is a multiset of first-order 00103 * derivatives, the order of the multiple derivative is the 00104 * size of the set. */ 00105 int order() const {return this->size();} 00106 00107 /** Return the order of spatial differentiation in this 00108 * derivative */ 00109 int spatialOrder() const ; 00110 00111 /** Return a multiindex representing the spatial derivs in this 00112 * multiple derivative */ 00113 MultiIndex spatialDeriv() const ; 00114 00115 00116 /** Return by reference argument the partitioning of 00117 * derivatives when this derivative is applied to a product. 00118 */ 00119 void productRulePermutations(ProductRulePerms& perms) const ; 00120 00121 /** return the product of two multiple derivatives, which 00122 * is the union of the two multisets */ 00123 MultipleDeriv product(const MultipleDeriv& other) const ; 00124 00125 00126 /** 00127 * Return a copy of this derivative, but with the specified single 00128 * derivative removed. For example, if <t>this</t> is 00129 * \f$ \{u,v,w\} \f$, calling <t>factorOutDeriv(v)</t> will 00130 * return \f$\{u,w\}\f$. 00131 */ 00132 MultipleDeriv factorOutDeriv(const Deriv& x) const ; 00133 00134 00135 /** 00136 * Return a copy of this derivative, but with all single derivs 00137 * in the the argument removed. For example, if <t>this</t> is 00138 * \f$ \{u,v,w\} \f$, calling <t>factorOutDeriv(u,v)</t> will 00139 * return \f$\{w\}\f$. 00140 */ 00141 MultipleDeriv factorOutDeriv(const MultipleDeriv& x) const ; 00142 00143 00144 /** 00145 * Return true is all single derivatives appearing in x also appear in this. 00146 */ 00147 bool containsDeriv(const MultipleDeriv& x) const ; 00148 00149 00150 /** 00151 * 00152 */ 00153 bool isInRequiredSet(const Set<MultiSet<int> >& funcCombinations, 00154 const Set<MultiIndex>& multiIndices) const ; 00155 00156 /** */ 00157 MultiSet<FunctionIdentifier> funcIDs() const ; 00158 00159 /** */ 00160 MultiSet<int> dofIDs() const ; 00161 00162 /** \name Utilities used in computing product rule permutations */ 00163 //@ 00164 /** Compute the n-th power of 2 */ 00165 static int pow2(int n); 00166 00167 /** Compute the n bits of an integer x < 2^n. The bits are ordered 00168 * least-to-most significant. 00169 */ 00170 static Array<int> bitsOfAnInteger(int x, int n); 00171 //@} 00172 private: 00173 }; 00174 00175 00176 /** 00177 * \brief Tranform the input set S as follows: for each multiple derivative 00178 * in S, apply the multiindex x to each of its component functional derivatives 00179 * in turn. Note that the multiindex may have negative indices. Results with 00180 * negative multiindices are ignored. 00181 */ 00182 Set<MultipleDeriv> applyTx(const Set<MultipleDeriv>& S, 00183 const MultiIndex& x); 00184 /** 00185 * \brief Filter the input set W, allowing only coordinate derivatives in the 00186 * direction x and functional derivatives whose associated functions 00187 * have nonzero evaluation points. This function is used 00188 * in the spatial/functional chain rule to identify those terms resulting 00189 * from differentiation wrt x or functional derivatives. 00190 */ 00191 Set<MultipleDeriv> applyZx(const Set<MultipleDeriv>& W, 00192 const MultiIndex& x); 00193 00194 00195 /** */ 00196 Set<MultipleDeriv> Xx(const MultiIndex& x) ; 00197 00198 /** */ 00199 int factorial(const MultipleDeriv& ms); 00200 00201 00202 /** 00203 * 00204 */ 00205 template <class T> class increasingOrder 00206 : std::binary_function<T, T, bool> 00207 { 00208 public: 00209 bool operator()(const MultipleDeriv& a, 00210 const MultipleDeriv& b) const; 00211 }; 00212 /** 00213 * Functor increasingOrder() is used as a comparison method 00214 * for MultipleDeriv objects. When used in a set, it will sort 00215 * derivatives in order of increasing order of differentiation. 00216 * Sorting by differentiation order is useful in evaluation of 00217 * product and chain rules: if we evaluate higher-order derivatives 00218 * first, we can evaluate in place without destroying lower-order 00219 * derivatives. 00220 */ 00221 template <> class increasingOrder<MultipleDeriv> 00222 { 00223 public: 00224 bool operator()(const MultipleDeriv& a, 00225 const MultipleDeriv& b) const 00226 { 00227 if (a.order() < b.order()) return true; 00228 if (a.order() > b.order()) return false; 00229 00230 MultipleDeriv::const_iterator aIter; 00231 MultipleDeriv::const_iterator bIter; 00232 00233 for (aIter=a.begin(), bIter=b.begin(); 00234 aIter != a.end() && bIter != b.end(); 00235 aIter++, bIter++) 00236 { 00237 if (*aIter < *bIter) return true; 00238 if (*bIter < *aIter) return false; 00239 } 00240 00241 return false; 00242 } 00243 }; 00244 00245 00246 /** */ 00247 MultipleDeriv makeMultiDeriv(const Deriv& d); 00248 00249 /** */ 00250 bool hasParameter(const MultipleDeriv& d) ; 00251 00252 00253 } 00254 00255 00256 #endif