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_EXPR_H 00032 #define SUNDANCE_EXPR_H 00033 00034 #include "SundanceDefs.hpp" 00035 #include "SundanceExprBase.hpp" 00036 #include "SundanceFunctionIdentifier.hpp" 00037 #include "SundanceMap.hpp" 00038 #include "PlayaHandle.hpp" 00039 #include "Teuchos_RefCountPtr.hpp" 00040 #include "Teuchos_TimeMonitor.hpp" 00041 #include <complex> 00042 00043 00044 namespace Sundance 00045 { 00046 00047 class ScalarExpr; 00048 class FuncElementBase; 00049 00050 using namespace Sundance; 00051 /** 00052 * User-level expression class. Expr is a handle to a 00053 * reference-counted pointer to a ExprBase subtype. As such, 00054 * expression copies and assignments are shallow. 00055 * 00056 * <h2> Lists </h2> 00057 * 00058 * Expressions can be grouped into lists with an arbitrary structure. Important 00059 * special cases are scalar, vector, and tensor expressions. 00060 * 00061 * <h3> Creating Lists </h3> 00062 * 00063 * <ul> 00064 * <li> Vector expressions 00065 * \code 00066 * Expr v = List(vx, vy, vz); 00067 * \endcode 00068 * <li> Heterogeneous lists 00069 * \code 00070 * // Form vector {vx, vy, vz} 00071 * Expr v = List(vx, vy, vz); 00072 * Expr q = new TestFunction(new Lagrange(1)); 00073 * // Form heterogeneous list {{vx, vy, vz}, q} 00074 * Expr state = List(v, q); 00075 * \endcode 00076 * </ul> 00077 * 00078 * <h3> Probing Lists </h3> 00079 * 00080 * <ul> 00081 * <li> Finding the size of the top level of a list 00082 * \code 00083 * // Form vector {vx, vy, vz} 00084 * Expr v = List(vx, vy, vz); 00085 * Expr q = new TestFunction(new Lagrange(1)); 00086 * // Form heterogeneous list {{vx, vy, vz}, q} 00087 * Expr state = List(v, q); 00088 * // Check top-level size of state list {{vx, vy, vz}, q} 00089 * int stateSize = state.size(); // returns 2 00090 * \endcode 00091 * <li> Finding the total size of a list 00092 * \code 00093 * // Check total size of state list 00094 * int totalSize = state.totalSize(); // returns 4 00095 * \endcode 00096 * </ul> 00097 * 00098 * <h3> Manipulating Lists </h3> 00099 * 00100 * 00101 * 00102 * <h2> Arithmetic Operations </h2> 00103 * 00104 * <ul> 00105 * <li> Addition 00106 * \code 00107 * Expr sum = x + y; 00108 * \endcode 00109 * The operands must have identical list structures. 00110 * 00111 * <li> Subtraction 00112 * \code 00113 * Expr diff = x - y; 00114 * \endcode 00115 * The operands must have identical list structures. 00116 * 00117 * <li> Multiplication 00118 * \code 00119 * Expr product = x * y; 00120 * \endcode 00121 * The operands must have list 00122 * structures such that the product can be interpreted as a 00123 * scalar-vector product or as an inner product between vectors 00124 * or tensors. The multiplication operator is also used to 00125 * represent the application of a differential operator. 00126 * 00127 * <li> Division 00128 * \code 00129 * Expr quotient = x / y; 00130 * \endcode 00131 * The denominator must be scalar-valued. 00132 * </ul> 00133 * 00134 * <h2> Expression Subtypes </h2> 00135 * The user-level expression subtypes are listed below, along with examples of their use. 00136 * <ul> 00137 * <li> UnknownFunction - Represents an unknown function in a finite-element problem. 00138 * Unknown functions can be scalar-valued or vector valued. 00139 * 00140 * Example of creation of a scalar-valued unknown function: 00141 * \code 00142 * Expr u = new UnknownFunction(new Lagrange(1)); 00143 * \endcode 00144 * 00145 * Example of creation of a vector-valued unknown function: 00146 * \code 00147 * Expr u = new UnknownFunction(List(new Lagrange(1), new Lagrange(1))); 00148 * \endcode 00149 * 00150 * <li> TestFunction - Represents a test function in a finite-element problem. 00151 * Test functions can be scalar-valued or vector valued. 00152 * 00153 * Example of creation of a scalar-valued test function: 00154 * \code 00155 * Expr v = new TestFunction(new Lagrange(1)); 00156 * \endcode 00157 * 00158 * Example of creation of a vector-valued test function: 00159 * \code 00160 * Expr u = new TestFunction(List(new Lagrange(1), new Lagrange(1))); 00161 * \endcode 00162 * 00163 * <li> Derivative - Represents a spatial derivative operator. 00164 * Spatial derivatives are applied 00165 * using the multiplication operator. 00166 * \code 00167 * Expr dx = new Derivative(0); 00168 * Expr convect = (u*dx)*u; 00169 * \endcode 00170 * Derivative expressions are scalar valued. However, vector differential operators 00171 * can be created using the List operator. For example, 00172 * \code 00173 * Expr dx = new Derivative(0); 00174 * Expr dy = new Derivative(1); 00175 * Expr grad = List(dx, dy); 00176 * \endcode 00177 * 00178 * <li> CoordExpr - Represents a coordinate functions. 00179 * \code 00180 * Expr x = new CoordExpr(0); 00181 * Expr y = new CoordExpr(1); 00182 * Expr r = sqrt(x*x + y*y); 00183 * \endcode 00184 * Coordinate expressions are scalar valued. 00185 * 00186 * <li> CellDiameterExpr - Represents the diameter of a cell. Cell 00187 * diameters are often used in stabilized methods. 00188 * \code 00189 * Expr h = new CellDiameterExpr(); 00190 * Expr streamlineDiffusion = eps*h*h*((u*grad)*v)*((u*grad)*T); 00191 * \endcode 00192 * Cell diameter expressions are scalar valued. 00193 * </ul> 00194 */ 00195 class Expr : public Playa::Handle<ExprBase> 00196 { 00197 public: 00198 /* boilerplate handle ctors */ 00199 HANDLE_CTORS(Expr, ExprBase); 00200 00201 /** Construct with a constant. Creates a ConstantExpr. */ 00202 Expr(const double& c); 00203 00204 /** Construct with a complex constant. Creates a ComplexExpr. */ 00205 Expr(const std::complex<double>& c); 00206 00207 /** Add two expressions. The operands must have identical list 00208 * structures. 00209 */ 00210 Expr operator+(const Expr& other) const ; 00211 /** Subtract two expressions. The operands must have identical 00212 * list structures. */ 00213 Expr operator-(const Expr& other) const ; 00214 /** Multiply two expressions. The operands must have list 00215 * structures such that the product can be interpreted as a 00216 * scalar-vector product or as an inner product between vectors 00217 * or tensors. The multiplication operator is also used to 00218 * represent the application of a differential operator. 00219 */ 00220 Expr operator*(const Expr& other) const ; 00221 /** Divide one expression by another. The right operand must be 00222 a scalar. */ 00223 Expr operator/(const Expr& other) const ; 00224 00225 /** Divide an expression by a double. */ 00226 Expr operator/(const double& other) const 00227 {return operator*(1.0/other);} 00228 00229 /** Divide an expression by a complex. */ 00230 Expr operator/(const std::complex<double>& other) const 00231 {return operator*(1.0/other);} 00232 00233 /** Unary minus operator */ 00234 Expr operator-() const ; 00235 00236 /** Return real part of a complex expression */ 00237 Expr real() const ; 00238 00239 /** Return imaginary part of a complex expression */ 00240 Expr imag() const ; 00241 00242 /** Return complex conjugate */ 00243 Expr conj() const ; 00244 00245 00246 /** List element accessor */ 00247 const Expr& operator[](int i) const ; 00248 00249 /** Number of elements in top level of list */ 00250 int size() const ; 00251 00252 /** Total number of elements in list. */ 00253 int totalSize() const ; 00254 00255 /** Append a new element to this list */ 00256 void append(const Expr& expr); 00257 00258 /** Flatten this list */ 00259 Expr flatten() const ; 00260 00261 /** Flatten list and spectral structure */ 00262 Expr flattenSpectral() const ; 00263 00264 /** */ 00265 std::string toString() const ; 00266 00267 /** */ 00268 XMLObject toXML() const ; 00269 00270 /** */ 00271 bool isComplex() const ; 00272 00273 /** */ 00274 bool isSpectral() const; 00275 00276 /** */ 00277 bool isTestElement() const ; 00278 00279 /** */ 00280 bool isUnknownElement() const ; 00281 00282 /** */ 00283 void setParameterValue(const double& value); 00284 /** */ 00285 double getParameterValue() const ; 00286 00287 /** Indicate whether the expression is independent of the given 00288 * functions */ 00289 bool isIndependentOf(const Expr& u) const ; 00290 00291 00292 /** 00293 * Indicate whether the expression is nonlinear 00294 * with respect to test functions */ 00295 bool isLinearInTests() const ; 00296 00297 /** 00298 * Indicate whether every term in the expression contains test functions */ 00299 bool everyTermHasTestFunctions() const ; 00300 00301 /** 00302 * Indicate whether the expression contains test functions */ 00303 bool hasTestFunctions() const ; 00304 00305 /** Indicate whether the expression is linear in the given 00306 * functions */ 00307 bool isLinearForm(const Expr& u) const ; 00308 00309 /** Indicate whether the expression is quadratic in the given 00310 * functions */ 00311 bool isQuadraticForm(const Expr& u) const ; 00312 00313 /** Find the unknown functions appearing in this expression */ 00314 void getUnknowns(Set<int>& unkID, Array<Expr>& unks) const ; 00315 00316 /** Find the test functions appearing in this expression */ 00317 void getTests(Set<int>& testID, Array<Expr>& tests) const ; 00318 00319 00320 #ifndef DOXYGEN_DEVELOPER_ONLY 00321 00322 00323 00324 /** 00325 * Turn evaluation caching on 00326 */ 00327 static bool& evaluationCachingOn() {static bool rtn = true; return rtn;} 00328 00329 /** 00330 * Show parentheses around every pair of operands 00331 */ 00332 static bool& showAllParens() {static bool rtn = false; return rtn;} 00333 00334 /** Create a new handle for an existing ptr. */ 00335 static Expr handle(const RCP<ExprBase>& ptr) 00336 {return Expr(ptr);} 00337 00338 /** */ 00339 static Time& opTimer() 00340 { 00341 static RCP<Time> rtn 00342 = TimeMonitor::getNewTimer("Expr symbolic ops"); 00343 return *rtn; 00344 } 00345 /** */ 00346 static Time& outputTimer() 00347 { 00348 static RCP<Time> rtn 00349 = TimeMonitor::getNewTimer("Expr output"); 00350 return *rtn; 00351 } 00352 00353 /** Comparison operator for ordering expr trees, DO NOT use for comparison 00354 * of numerical values. */ 00355 bool lessThan(const Expr& other) const ; 00356 00357 /** Comparison operator for ordering expr trees, DO NOT use for comparison 00358 * of numerical values. */ 00359 bool operator<(const Expr& other) const ; 00360 00361 /** Equality test operator for ordering expr trees, DO NOT use for comparison 00362 * of numerical values. */ 00363 bool sameAs(const Expr& other) const ; 00364 00365 /** */ 00366 Sundance::Map<Expr, int> getSumTree() const ; 00367 00368 00369 00370 /** */ 00371 const ScalarExpr* scalarExpr() const ; 00372 00373 /** */ 00374 const FuncElementBase* funcElement() const ; 00375 00376 private: 00377 00378 00379 00380 /** Add two scalar expressions */ 00381 Expr sum(const Expr& other, int sign) const ; 00382 00383 /** Multiply two scalar expressions */ 00384 Expr multiply(const Expr& other) const ; 00385 00386 /** Divide two scalar expressions */ 00387 Expr divide(const Expr& other) const ; 00388 00389 /** Try transformations of complex addition */ 00390 bool tryAddComplex(const Expr& L, const Expr& R, int sign, 00391 Expr& rtn) const ; 00392 00393 /** Try transformations of complex multiplication */ 00394 bool tryMultiplyComplex(const Expr& L, const Expr& R, 00395 Expr& rtn) const ; 00396 00397 00398 00399 00400 #endif /* DOXYGEN_DEVELOPER_ONLY */ 00401 }; 00402 00403 /** \relates Expr */ 00404 inline std::ostream& operator<<(std::ostream& os, const Expr& e) 00405 { 00406 if (e.ptr().get()==0) {os << "Expr()"; return os;} 00407 return e.ptr()->toText(os, false); 00408 } 00409 00410 /** \relates Expr */ 00411 inline Expr operator+(const double& a, const Expr& x) 00412 {return Expr(a) + x;} 00413 00414 /** \relates Expr */ 00415 inline Expr operator-(const double& a, const Expr& x) 00416 {return Expr(a) - x;} 00417 00418 /** \relates Expr */ 00419 inline Expr operator*(const double& a, const Expr& x) 00420 {return Expr(a) * x;} 00421 00422 /** \relates Expr */ 00423 inline Expr operator/(const double& a, const Expr& x) 00424 {return Expr(a) / x;} 00425 00426 /** \relates Expr */ 00427 inline Expr operator+(const std::complex<double>& a, const Expr& x) 00428 {return Expr(a) + x;} 00429 00430 /** \relates Expr */ 00431 inline Expr operator-(const std::complex<double>& a, const Expr& x) 00432 {return Expr(a) - x;} 00433 00434 /** \relates Expr */ 00435 inline Expr operator*(const std::complex<double>& a, const Expr& x) 00436 {return Expr(a) * x;} 00437 00438 /** \relates Expr */ 00439 inline Expr operator/(const std::complex<double>& a, const Expr& x) 00440 {return Expr(a) / x;} 00441 00442 /** \relates Expr */ 00443 inline Expr Re(const Expr& a) {return a.real();} 00444 00445 /** \relates Expr */ 00446 inline Expr Im(const Expr& a) {return a.imag();} 00447 00448 /** \relates Expr */ 00449 inline Expr conj(const Expr& a) {return a.conj();} 00450 00451 /** \relates Expr */ 00452 Expr Complex(const Expr& real, const Expr& imag); 00453 00454 /** \relates Expr */ 00455 Expr List(const Expr& a); 00456 00457 /** \relates Expr */ 00458 Expr List(const Expr& a, const Expr& b); 00459 00460 /** \relates Expr */ 00461 Expr List(const Expr& a, const Expr& b, const Expr& c); 00462 00463 /** \relates Expr */ 00464 Expr List(const Expr& a, const Expr& b, const Expr& c, 00465 const Expr& d); 00466 00467 /** \relates Expr */ 00468 Expr List(const Expr& a, const Expr& b, const Expr& c, 00469 const Expr& d, const Expr& e); 00470 00471 /** \relates Expr */ 00472 Expr List(const Expr& a, const Expr& b, const Expr& c, 00473 const Expr& d, const Expr& e, const Expr& f); 00474 00475 /** \relates Expr */ 00476 Expr List(const Expr& a, const Expr& b, const Expr& c, 00477 const Expr& d, const Expr& e, const Expr& f, 00478 const Expr& g); 00479 00480 00481 /** \relates Expr */ 00482 Expr List(const Expr& a, const Expr& b, const Expr& c, 00483 const Expr& d, const Expr& e, const Expr& f, 00484 const Expr& g, const Expr& h); 00485 00486 /** \relates Expr */ 00487 Expr List(const Expr& a, const Expr& b, const Expr& c, 00488 const Expr& d, const Expr& e, const Expr& f, 00489 const Expr& g, const Expr& h, const Expr& i); 00490 00491 /** \relates Expr */ 00492 Expr List(const Expr& a, const Expr& b, const Expr& c, 00493 const Expr& d, const Expr& e, const Expr& f, 00494 const Expr& g, const Expr& h, const Expr& i, 00495 const Expr& j); 00496 00497 00498 00499 00500 } 00501 00502 #endif