SyFi 0.3
|
00001 // Copyright (C) 2006-2009 Kent-Andre Mardal and Simula Research Laboratory. 00002 // Licensed under the GNU GPL Version 2, or (at your option) any later version. 00003 00004 #include "ginac_tools.h" 00005 #include "symbol_factory.h" 00006 00007 #include <ginac/ginac.h> 00008 00009 #include <iostream> 00010 #include <sstream> 00011 #include <fstream> 00012 #include <stdexcept> 00013 00014 using namespace std; 00015 using namespace GiNaC; 00016 00017 namespace SyFi 00018 { 00019 00020 GiNaC::lst cross(GiNaC::lst& v1, GiNaC::lst& v2) 00021 { 00022 GiNaC::lst ret; 00023 if ( v1.nops() != v2.nops() ) 00024 { 00025 cout <<"incompatible vectors "<<endl; 00026 cout <<"v1.nops() "<<v1.nops(); 00027 cout <<" v2.nops() "<<v2.nops()<<endl; ; 00028 return GiNaC::lst(); 00029 } 00030 ret.append( v1.op(1)*v2.op(2) - v1.op(2)*v2.op(1)); 00031 ret.append(- v1.op(0)*v2.op(2) + v1.op(2)*v2.op(0)); 00032 ret.append( v1.op(0)*v2.op(1) - v1.op(1)*v2.op(0)); 00033 return ret; 00034 } 00035 00036 GiNaC::ex inner(GiNaC::ex a, GiNaC::ex b, bool transposed) 00037 { 00038 if (GiNaC::is_a<GiNaC::matrix>(a) && GiNaC::is_a<GiNaC::matrix>(b)) 00039 { 00040 GiNaC::matrix ma = GiNaC::ex_to<GiNaC::matrix>(a); 00041 GiNaC::matrix mb = GiNaC::ex_to<GiNaC::matrix>(b); 00042 if ( !transposed ) 00043 { 00044 if (ma.cols() != mb.cols() || ma.rows() != mb.rows() ) 00045 { 00046 cout <<"Incompatible matrices "<<endl; 00047 cout <<"a.cols() "<<ma.cols()<<endl; 00048 cout <<"a.rows() "<<ma.rows()<<endl; 00049 cout <<"b.cols() "<<mb.cols()<<endl; 00050 cout <<"b.rows() "<<mb.rows()<<endl; 00051 cout <<"a="<<a<<endl; 00052 cout <<"b="<<b<<endl; 00053 throw std::runtime_error("Incompatible matrices."); 00054 } 00055 00056 GiNaC::ex ret; 00057 for (unsigned int i=0; i<ma.rows(); i++) 00058 { 00059 for (unsigned int j=0; j<ma.cols(); j++) 00060 { 00061 ret += ma(i,j)*mb(i,j); 00062 } 00063 } 00064 return ret; 00065 } 00066 else 00067 { 00068 if (ma.cols() != mb.rows() || ma.rows() != mb.cols() ) 00069 { 00070 cout <<"Incompatible matrices "<<endl; 00071 cout <<"a.cols() "<<ma.cols()<<endl; 00072 cout <<"a.rows() "<<ma.rows()<<endl; 00073 cout <<"b.cols() "<<mb.cols()<<endl; 00074 cout <<"b.rows() "<<mb.rows()<<endl; 00075 cout <<"a="<<a<<endl; 00076 cout <<"b="<<b<<endl; 00077 throw std::runtime_error("Incompatible matrices."); 00078 } 00079 00080 GiNaC::ex ret; 00081 for (unsigned int i=0; i<ma.rows(); i++) 00082 { 00083 for (unsigned int j=0; j<ma.cols(); j++) 00084 { 00085 ret += ma(i,j)*mb(j,i); 00086 } 00087 } 00088 return ret; 00089 } 00090 } 00091 else if (GiNaC::is_a<GiNaC::lst>(a) 00092 && GiNaC::is_a<GiNaC::lst>(b)) 00093 { 00094 return inner(GiNaC::ex_to<GiNaC::lst>(a), GiNaC::ex_to<GiNaC::lst>(b)); 00095 } 00096 else 00097 { 00098 return a*b; 00099 } 00100 } 00101 00102 GiNaC::ex inner(GiNaC::lst v1, GiNaC::lst v2) 00103 { 00104 GiNaC::ex ret; 00105 00106 if ( v1.nops() != v2.nops() ) 00107 { 00108 cout <<"incompatible vectors "<<endl; 00109 cout <<"v1.nops() "<<v1.nops(); 00110 cout <<" v2.nops() "<<v2.nops()<<endl; ; 00111 return 0; 00112 } 00113 for (unsigned i = 0; i <= v1.nops()-1 ; ++i) 00114 { 00115 if ( GiNaC::is_a<GiNaC::lst>(v1.op(i)) && 00116 GiNaC::is_a<GiNaC::lst>(v2.op(i)) ) 00117 { 00118 ret += inner(GiNaC::ex_to<GiNaC::lst>(v1.op(i)), 00119 GiNaC::ex_to<GiNaC::lst>(v2.op(i))); 00120 } 00121 else 00122 { 00123 ret += v1.op(i)*v2.op(i); 00124 } 00125 } 00126 return ret; 00127 } 00128 00129 GiNaC::ex inner(GiNaC::exvector& v1, GiNaC::exvector& v2) 00130 { 00131 GiNaC::ex ret; 00132 for (unsigned int i=0; i< v1.size(); i++) 00133 { 00134 ret += v1[i]*v2[i]; 00135 } 00136 return ret; 00137 } 00138 00139 GiNaC::lst matvec(GiNaC::matrix& M, GiNaC::lst& x) 00140 { 00141 GiNaC::lst ret; 00142 int nr = M.rows(); 00143 int nc = M.cols(); 00144 for (int i = 0; i < nr; i++) 00145 { 00146 GiNaC::ex tmp; 00147 for (int j = 0; j < nc; j++) 00148 { 00149 tmp = tmp + M(i,j)*(x.op(j)); 00150 } 00151 ret.append(tmp); 00152 } 00153 return ret; 00154 } 00155 00156 GiNaC::ex matvec(GiNaC::ex A, GiNaC::ex x) 00157 { 00158 ex sol; 00159 00160 if (GiNaC::is_a<GiNaC::matrix>(A) && GiNaC::is_a<GiNaC::matrix>(x)) 00161 { 00162 GiNaC::matrix AA = GiNaC::ex_to<GiNaC::matrix>(A); 00163 GiNaC::matrix xx = GiNaC::ex_to<GiNaC::matrix>(x); 00164 sol = AA.mul(xx); 00165 } 00166 else 00167 { 00168 throw std::runtime_error("Invalid argument types, need matrices"); 00169 } 00170 return sol; 00171 } 00172 00173 GiNaC::lst ex2equations(GiNaC::ex rel) 00174 { 00175 GiNaC::ex lhs = rel.lhs(); 00176 GiNaC::ex rhs = rel.rhs(); 00177 00178 GiNaC::ex l; 00179 GiNaC::ex r; 00180 00181 GiNaC::lst eqs; 00182 00183 for (int i=lhs.ldegree(x); i<=lhs.degree(x); ++i) 00184 { 00185 for (int j=lhs.ldegree(y); j<=lhs.degree(y); ++j) 00186 { 00187 for (int k=lhs.ldegree(z); k<=lhs.degree(z); ++k) 00188 { 00189 l = lhs.coeff(x,i).coeff(y, j).coeff(z,k); 00190 r = rhs.coeff(x,i).coeff(y, j).coeff(z,k); 00191 // if (! (l == 0 && r == 0 ) ) eqs.append(l == r); OLD VERSION 00192 if ( (l != 0 && (r == 0 || r == 1) ) ) eqs.append(l == r); 00193 } 00194 } 00195 } 00196 eqs.sort(); 00197 return eqs; 00198 } 00199 00200 GiNaC::lst collapse(GiNaC::lst l) 00201 { 00202 GiNaC::lst lc; 00203 GiNaC::lst::const_iterator iter1, iter2; 00204 00205 for (iter1 = l.begin(); iter1 != l.end(); ++iter1) 00206 { 00207 if (GiNaC::is_a<GiNaC::lst>(*iter1)) 00208 { 00209 for (iter2 = GiNaC::ex_to<GiNaC::lst>(*iter1).begin(); iter2 != GiNaC::ex_to<GiNaC::lst>(*iter1).end(); ++iter2) 00210 { 00211 lc.append(*iter2); 00212 } 00213 } 00214 else 00215 { 00216 lc.append(*iter1); 00217 } 00218 } 00219 lc.sort(); 00220 lc.unique(); 00221 return lc; 00222 } 00223 00224 GiNaC::matrix equations2matrix(const GiNaC::ex &eqns, const GiNaC::ex &symbols) 00225 { 00226 00227 GiNaC::matrix sys(eqns.nops(),symbols.nops()); 00228 GiNaC::matrix rhs(eqns.nops(),1); 00229 GiNaC::matrix vars(symbols.nops(),1); 00230 00231 for (size_t r=0; r<eqns.nops(); r++) 00232 { 00233 // lhs-rhs==0 00234 const GiNaC::ex eq = eqns.op(r).op(0)-eqns.op(r).op(1); 00235 GiNaC::ex linpart = eq; 00236 for (size_t c=0; c<symbols.nops(); c++) 00237 { 00238 const GiNaC::ex co = eq.coeff(GiNaC::ex_to<GiNaC::symbol>(symbols.op(c)),1); 00239 linpart -= co*symbols.op(c); 00240 sys(r,c) = co; 00241 } 00242 linpart = linpart.expand(); 00243 rhs(r,0) = -linpart; 00244 } 00245 return sys; 00246 } 00247 00248 void matrix_from_equations(const GiNaC::ex &eqns, const GiNaC::ex &symbols, GiNaC::matrix& A, GiNaC::matrix& b) 00249 { 00250 // build matrix from equation system 00251 GiNaC::matrix sys(eqns.nops(),symbols.nops()); 00252 GiNaC::matrix rhs(eqns.nops(),1); 00253 GiNaC::matrix vars(symbols.nops(),1); 00254 00255 for (size_t r=0; r<eqns.nops(); r++) 00256 { 00257 // lhs-rhs==0 00258 const GiNaC::ex eq = eqns.op(r).op(0)-eqns.op(r).op(1); 00259 GiNaC::ex linpart = eq; 00260 for (size_t c=0; c<symbols.nops(); c++) 00261 { 00262 const GiNaC::ex co = eq.coeff(GiNaC::ex_to<GiNaC::symbol>(symbols.op(c)),1); 00263 linpart -= co*symbols.op(c); 00264 sys(r,c) = co; 00265 } 00266 linpart = linpart.expand(); 00267 rhs(r,0) = -linpart; 00268 } 00269 A = sys; 00270 b = rhs; 00271 } 00272 00273 GiNaC::ex lst_to_matrix2(const GiNaC::lst& l) 00274 { 00275 GiNaC::lst::const_iterator itr, itc; 00276 00277 // Find number of rows and columns 00278 size_t rows = l.nops(), cols = 0; 00279 for (itr = l.begin(); itr != l.end(); ++itr) 00280 { 00281 if (!GiNaC::is_a<GiNaC::lst>(*itr)) 00282 // throw (std::invalid_argument("lst_to_matrix: argument must be a list of lists")); 00283 cols = 1; 00284 if (itr->nops() > cols) 00285 cols = itr->nops(); 00286 } 00287 // Allocate and fill matrix 00288 GiNaC::matrix &M = *new GiNaC::matrix(rows, cols); 00289 M.setflag(GiNaC::status_flags::dynallocated); 00290 00291 unsigned i; 00292 for (itr = l.begin(), i = 0; itr != l.end(); ++itr, ++i) 00293 { 00294 unsigned j; 00295 if (cols == 1) 00296 { 00297 M(i, 0) = *itr; 00298 } 00299 else 00300 { 00301 for (itc = GiNaC::ex_to<GiNaC::lst>(*itr).begin(), j = 0; itc != GiNaC::ex_to<GiNaC::lst>(*itr).end(); ++itc, ++j) 00302 M(i, j) = *itc; 00303 } 00304 } 00305 return M; 00306 } 00307 00308 GiNaC::lst matrix_to_lst2(const GiNaC::ex& m) 00309 { 00310 if (GiNaC::is_a<GiNaC::matrix>(m)) 00311 { 00312 GiNaC::matrix A = GiNaC::ex_to<GiNaC::matrix>(m); 00313 int cols = A.cols(); 00314 int rows = A.rows(); 00315 00316 GiNaC::lst ret; 00317 if ( cols == 1 ) 00318 { 00319 for (unsigned int i=0; i<=A.rows()-1; i++) 00320 { 00321 ret.append(A(i,0)); 00322 } 00323 } 00324 else if ( rows == 1 ) 00325 { 00326 for (unsigned int i=0; i<=A.cols()-1; i++) 00327 { 00328 ret.append(A(0,i)); 00329 } 00330 } 00331 else 00332 { 00333 for (unsigned int i=0; i<=A.rows()-1; i++) 00334 { 00335 GiNaC::lst rl; 00336 for (unsigned int j=0; j<=A.cols()-1; j++) 00337 { 00338 rl.append(A(i,j)); 00339 } 00340 ret.append(rl); 00341 } 00342 } 00343 return ret; 00344 } 00345 else 00346 { 00347 return GiNaC::lst(); 00348 } 00349 } 00350 00351 GiNaC::lst lst_equals(GiNaC::ex a, GiNaC::ex b) 00352 { 00353 GiNaC::lst ret; 00354 if ( (GiNaC::is_a<GiNaC::lst>(a)) && (GiNaC::is_a<GiNaC::lst>(b)) /*&& (a.nops() == b.nops())*/ ) { 00355 for (unsigned int i=0; i<= a.nops()-1; i++) 00356 { 00357 ret.append(b.op(i) == a.op(i)); 00358 } 00359 } 00360 else if ( !(GiNaC::is_a<GiNaC::lst>(a)) && !(GiNaC::is_a<GiNaC::lst>(b))) 00361 { 00362 ret.append(b == a); 00363 } 00364 else if ( !(GiNaC::is_a<GiNaC::lst>(a)) && (GiNaC::is_a<GiNaC::lst>(b))) 00365 { 00366 ret.append(b.op(0) == a); 00367 } 00368 else 00369 { 00370 throw(std::invalid_argument("Make sure that the lists a and b are comparable.")); 00371 } 00372 return ret; 00373 } 00374 00375 00376 int find(GiNaC::ex e, GiNaC::lst list) 00377 { 00378 for (unsigned int i=0; i< list.nops(); i++) 00379 { 00380 if ( e == list.op(i) ) return i; 00381 } 00382 return -1; 00383 } 00384 00385 00386 void visitor_subst_pow(GiNaC::ex e, GiNaC::exmap& map, ex_int_map& intmap, string a) 00387 { 00388 static int i=0; 00389 if (map.find(e) != map.end()) 00390 { 00391 intmap[e] = intmap[e]+1; 00392 return; 00393 } 00394 if (GiNaC::is_exactly_a<GiNaC::power>(e)) 00395 { 00396 std::ostringstream s; 00397 s <<a<<i++; 00398 map[e] = GiNaC::symbol(s.str()); 00399 intmap[e] = 0; 00400 for (unsigned int i=0; i< e.nops(); i++) 00401 { 00402 GiNaC::ex e2 = e.op(i); 00403 // cout <<"power e "<<e2<<endl; 00404 visitor_subst_pow(e2,map,intmap, a); 00405 } 00406 } 00407 else if (GiNaC::is_a<GiNaC::function>(e)) 00408 { 00409 std::ostringstream s; 00410 s <<a<<i++; 00411 map[e] = GiNaC::symbol(s.str()); 00412 intmap[e] = 0; 00413 for (unsigned int i=0; i< e.nops(); i++) 00414 { 00415 GiNaC::ex e2 = e.op(i); 00416 // cout <<"function e "<<e2<<endl; 00417 visitor_subst_pow(e2,map,intmap, a); 00418 } 00419 } 00420 else if (GiNaC::is_a<GiNaC::mul>(e)) 00421 { 00422 if (e.nops() > 4 && e.nops() < 10 ) 00423 { 00424 std::ostringstream s; 00425 s <<a<<i++; 00426 map[e] = GiNaC::symbol(s.str()); 00427 intmap[e] = 0; 00428 } 00429 00430 for (unsigned int i=0; i< e.nops(); i++) 00431 { 00432 GiNaC::ex e2 = e.op(i); 00433 visitor_subst_pow(e2,map,intmap, a); 00434 } 00435 } 00436 else if (GiNaC::is_a<GiNaC::add>(e)) 00437 { 00438 for (unsigned int i=0; i< e.nops(); i++) 00439 { 00440 GiNaC::ex e2 = e.op(i); 00441 visitor_subst_pow(e2,map,intmap,a); 00442 } 00443 } 00444 } 00445 00446 00447 void check_visitor(GiNaC::ex e, GiNaC::lst& exlist) 00448 { 00449 if (find(e, exlist) >= 0) return; 00450 00451 // cout <<"ex e "<<e<<endl; 00452 if (GiNaC::is_a<GiNaC::numeric>(e)) 00453 { 00454 } 00455 else if (GiNaC::is_a<GiNaC::add>(e) ) 00456 { 00457 // cout <<"e "<<e <<endl; 00458 // cout <<"e.nops() "<<e.nops() <<endl; 00459 if (e.nops() > 4 && e.nops() < 10 ) exlist.append(e); 00460 for (unsigned int i=0; i< e.nops(); i++) 00461 { 00462 GiNaC::ex e2 = e.op(i); 00463 // cout <<"add e "<<e2<<endl; 00464 // exlist.append(e2); 00465 check_visitor(e2,exlist); 00466 } 00467 } 00468 else if (GiNaC::is_a<GiNaC::mul>(e)) 00469 { 00470 for (unsigned int i=0; i< e.nops(); i++) 00471 { 00472 GiNaC::ex e2 = e.op(i); 00473 // cout <<"mul e "<<e2<<endl; 00474 exlist.append(e2); 00475 check_visitor(e2,exlist); 00476 } 00477 } 00478 else if (GiNaC::is_a<GiNaC::lst>(e)) 00479 { 00480 for (unsigned int i=0; i< e.nops(); i++) 00481 { 00482 GiNaC::ex e2 = e.op(i); 00483 // cout <<"GiNaC::lst e "<<e2<<endl; 00484 // exlist.append(e2); 00485 check_visitor(e2,exlist); 00486 } 00487 } 00488 else if (GiNaC::is_exactly_a<GiNaC::power>(e)) 00489 { 00490 exlist.append(e); 00491 for (unsigned int i=0; i< e.nops(); i++) 00492 { 00493 GiNaC::ex e2 = e.op(i); 00494 // cout <<"power e "<<e2<<endl; 00495 check_visitor(e2,exlist); 00496 } 00497 } 00498 else if (GiNaC::is_a<GiNaC::function>(e)) 00499 { 00500 exlist.append(e); 00501 for (unsigned int i=0; i< e.nops(); i++) 00502 { 00503 GiNaC::ex e2 = e.op(i); 00504 // cout <<"function e "<<e2<<endl; 00505 check_visitor(e2,exlist); 00506 } 00507 } 00508 00509 else 00510 { 00511 // exlist.append(e); 00512 // cout <<"atom e "<<e<<endl; 00513 } 00514 00515 exlist.sort(); 00516 exlist.unique(); 00517 } 00518 00519 00520 GiNaC::ex homogenous_pol(unsigned int order, unsigned int nsd, const string a) 00521 { 00522 using SyFi::x; 00523 using SyFi::y; 00524 using SyFi::z; 00525 00526 if ( nsd == 1) 00527 { 00528 GiNaC::symbol a0(istr(a,0)); 00529 return GiNaC::lst(a0*pow(x,order), a0, pow(x,order)); 00530 } 00531 else if ( nsd == 2 ) 00532 { 00533 GiNaC::ex variables = get_symbolic_matrix(1,order+1, a); 00534 GiNaC::lst basis; 00535 GiNaC::ex ret; 00536 for (unsigned int i=0; i<= order; i++) 00537 { 00538 basis.append(pow(x,i)*pow(y,order-i)); 00539 ret += variables.op(i)*basis.op(i); 00540 } 00541 return GiNaC::lst(ret, matrix_to_lst2(variables), basis); 00542 } 00543 else if ( nsd == 3 ) 00544 { 00545 GiNaC::lst basis; 00546 for (unsigned int i=0; i<= order; i++) 00547 { 00548 for (unsigned int j=0; j<= order; j++) 00549 { 00550 for (unsigned int k=0; k<= order; k++) 00551 { 00552 if ( i + j + k == order ) 00553 { 00554 basis.append(pow(x,i)*pow(y,j)*pow(z,k)); 00555 } 00556 } 00557 } 00558 } 00559 GiNaC::ex variables = get_symbolic_matrix(1,basis.nops(), a); 00560 GiNaC::ex ret; 00561 for (unsigned int i=0; i<basis.nops(); i++) 00562 { 00563 ret += variables.op(i)*basis.op(i); 00564 } 00565 return GiNaC::lst(ret, matrix_to_lst2(variables), basis); 00566 } 00567 throw std::runtime_error("Homogenous polynomials only implemented in 1D, 2D and 3D"); 00568 } 00569 00570 00571 GiNaC::lst homogenous_polv(unsigned int no_fields, unsigned int order, unsigned int nsd, const string a) 00572 { 00573 GiNaC::lst ret1; // contains the polynom 00574 GiNaC::lst ret2; // contains the coefficients 00575 GiNaC::lst ret3; // constains the basis functions 00576 GiNaC::lst basis_tmp; 00577 for (unsigned int i=0; i< no_fields; i++) 00578 { 00579 GiNaC::lst basis; 00580 std::ostringstream s; 00581 s <<a<<""<<i<<"_"; 00582 GiNaC::ex polspace = homogenous_pol(order, nsd, s.str()); 00583 ret1.append(polspace.op(0)); 00584 ret2.append(polspace.op(1)); 00585 basis_tmp = GiNaC::ex_to<GiNaC::lst>(polspace.op(2)); 00586 for (GiNaC::lst::const_iterator basis_iterator = basis_tmp.begin(); 00587 basis_iterator != basis_tmp.end(); ++basis_iterator) 00588 { 00589 GiNaC::lst tmp_lst; 00590 for (unsigned int d=1; d<=no_fields; d++) tmp_lst.append(0); 00591 tmp_lst.let_op(i) = (*basis_iterator); 00592 ret3.append(tmp_lst); 00593 } 00594 } 00595 return GiNaC::lst(ret1,ret2,ret3); 00596 } 00597 00598 00599 GiNaC::ex pol(unsigned int order, unsigned int nsd, const string a) 00600 { 00601 using SyFi::x; 00602 using SyFi::y; 00603 using SyFi::z; 00604 00605 GiNaC::ex ret; // ex to return 00606 int dof; // degrees of freedom 00607 GiNaC::ex A; // ex holding the coefficients a_0 .. a_dof 00608 GiNaC::lst basis; 00609 00610 if (nsd == 1) 00611 { 00612 /* 1D: 00613 * P^n = a_0 + a_1*x + .... + a_n*x^n 00614 * dof : n+1 00615 */ 00616 dof = order+1; 00617 A = get_symbolic_matrix(1,dof, a); 00618 int o=0; 00619 for (GiNaC::const_iterator i = A.begin(); i != A.end(); ++i) 00620 { 00621 ret += (*i)*pow(x,o); 00622 basis.append(pow(x,o)); 00623 o++; 00624 } 00625 } 00626 else if ( nsd == 2) 00627 { 00628 00629 /* 2D: structure of coefficients (a_i) 00630 * [ a_0 a_1 x a_3 x^2 a_6 x^3 00631 * [ a_2 y a_4 xy a_7 x^2y 00632 * [ a_5 y^2 a_8 xy^2 00633 * [ a_9 y^3 00634 */ 00635 dof = (order+1)*(order+2)/2; 00636 A = get_symbolic_matrix(1, dof , a); 00637 00638 size_t i=0; 00639 for (unsigned int o = 0; o <= order; o++) 00640 { 00641 for (unsigned int d = 0; d <= o; d++) 00642 { 00643 ret += A.op(i)*pow(y,d)*pow(x,o-d); 00644 basis.append(pow(y,d)*pow(x,o-d)); 00645 i++; 00646 } 00647 } 00648 } 00649 else if (nsd == 3) 00650 { 00651 00652 /* Similar structure as in 2D, but 00653 * structured as a tetraheder, i.e., 00654 * a_o + a_1 x + a_2 y + a_3 z 00655 * + a_4 x^2 + a_5 xy + 00656 */ 00657 dof = 0; 00658 for (unsigned int j=0; j<= order; j++) 00659 { 00660 dof += (j+1)*(j+2)/2; 00661 } 00662 A = get_symbolic_matrix(1, dof , a); 00663 00664 size_t i=0; 00665 for (unsigned int o = 0; o <= order; o++) 00666 { 00667 for (unsigned int d = 0; d <= o; d++) 00668 { 00669 for (unsigned int f = 0; f <= o; f++) 00670 { 00671 if ( int(o)-int(d)-int(f) >= 0) 00672 { 00673 ret += A.op(i)*pow(y,f)*pow(z,d)*pow(x,o-d-f); 00674 basis.append(pow(y,f)*pow(z,d)*pow(x,o-d-f)); 00675 i++; 00676 } 00677 } 00678 } 00679 } 00680 } 00681 return GiNaC::lst(ret,matrix_to_lst2(A), basis); 00682 } 00683 00684 00685 GiNaC::lst polv(unsigned int no_fields, unsigned int order, unsigned int nsd, const string a) 00686 { 00687 GiNaC::lst ret1; // contains the polynom 00688 GiNaC::lst ret2; // contains the coefficients 00689 GiNaC::lst ret3; // constains the basis functions 00690 GiNaC::lst basis_tmp; 00691 for (unsigned int i=0; i< no_fields; i++) 00692 { 00693 GiNaC::lst basis; 00694 std::ostringstream s; 00695 s <<a<<""<<i<<"_"; 00696 GiNaC::ex polspace = pol(order, nsd, s.str()); 00697 ret1.append(polspace.op(0)); 00698 ret2.append(polspace.op(1)); 00699 basis_tmp = GiNaC::ex_to<GiNaC::lst>(polspace.op(2)); 00700 for (GiNaC::lst::const_iterator basis_iterator = basis_tmp.begin(); 00701 basis_iterator != basis_tmp.end(); ++basis_iterator) 00702 { 00703 GiNaC::lst tmp_lst; 00704 for (unsigned int d=1; d<=no_fields; d++) tmp_lst.append(0); 00705 tmp_lst.let_op(i) = (*basis_iterator); 00706 ret3.append(tmp_lst); 00707 } 00708 } 00709 return GiNaC::lst(ret1,ret2,ret3); 00710 00711 /* Old Code: 00712 GiNaC::lst ret; 00713 for (int i=1; i<= nsd; i++) { 00714 std::ostringstream s; 00715 s <<a<<"^"<<i<<"_"; 00716 GiNaC::ex p = pol(order, nsd, s.str()); 00717 ret.append(p); 00718 } 00719 return ret; 00720 */ 00721 } 00722 00723 00724 GiNaC::ex polb(unsigned int order, unsigned int nsd, const string a) 00725 { 00726 using SyFi::x; 00727 using SyFi::y; 00728 using SyFi::z; 00729 00730 GiNaC::ex ret; // ex to return 00731 int dof; // degrees of freedom 00732 GiNaC::ex A; // ex holding the coefficients a_0 .. a_dof 00733 GiNaC::lst basis; 00734 00735 if (nsd == 1) 00736 { 00737 /* 1D: 00738 * P^n = a_0 + a_1*x + .... + a_n*x^n 00739 * dof : n+1 00740 */ 00741 dof = order+1; 00742 A = get_symbolic_matrix(1,dof, a); 00743 int o=0; 00744 for (GiNaC::const_iterator i = A.begin(); i != A.end(); ++i) 00745 { 00746 ret += (*i)*pow(x,o); 00747 basis.append(pow(x,o)); 00748 o++; 00749 } 00750 } 00751 else if ( nsd == 2) 00752 { 00753 00754 /* 2D: structure of coefficients (a_i) 00755 * [ a_0 a_1 x a_3 x^2 a_6 x^3 00756 * [ a_2 y a_4 xy a_7 x^2y 00757 * [ a_5 y^2 a_8 xy^2 00758 * [ a_9 y^3 00759 */ 00760 00761 dof = (order+1)*(order+1); 00762 A = get_symbolic_matrix(1, dof , a); 00763 00764 size_t i=0; 00765 for (unsigned int o = 0; o <= order; o++) 00766 { 00767 for (unsigned int d = 0; d <= order; d++) 00768 { 00769 ret += A.op(i)*pow(y,d)*pow(x,o); 00770 basis.append(pow(y,d)*pow(x,o)); 00771 i++; 00772 } 00773 } 00774 } 00775 else if (nsd == 3) 00776 { 00777 00778 /* Similar structure as in 2D, but 00779 * structured as a tetraheder, i.e., 00780 * a_o + a_1 x + a_2 y + a_3 z 00781 * + a_4 x^2 + a_5 xy + 00782 */ 00783 dof = (order+1)*(order+1)*(order+1); 00784 A = get_symbolic_matrix(1, dof , a); 00785 00786 size_t i=0; 00787 for (unsigned int o = 0; o <= order; o++) 00788 { 00789 for (unsigned int d = 0; d <= order; d++) 00790 { 00791 for (unsigned int f = 0; f <= order; f++) 00792 { 00793 ret += A.op(i)*pow(y,f)*pow(z,d)*pow(x,o); 00794 basis.append(pow(y,f)*pow(z,d)*pow(x,o)); 00795 i++; 00796 } 00797 } 00798 } 00799 } 00800 00801 return GiNaC::lst(ret,matrix_to_lst2(A), basis); 00802 } 00803 00804 00805 GiNaC::lst coeffs(GiNaC::lst pols) 00806 { 00807 GiNaC::lst cc; 00808 GiNaC::lst tmp; 00809 for (unsigned int i=0; i<= pols.nops()-1; i++) 00810 { 00811 tmp = coeffs(pols.op(i)); 00812 cc = collapse(GiNaC::lst(cc, tmp)); 00813 } 00814 return cc; 00815 } 00816 00817 00818 GiNaC::lst coeffs(GiNaC::ex pol) 00819 { 00820 using SyFi::x; 00821 using SyFi::y; 00822 using SyFi::z; 00823 00824 GiNaC::lst cc; 00825 GiNaC::ex c, b; 00826 for (int i=pol.ldegree(x); i<=pol.degree(x); ++i) 00827 { 00828 for (int j=pol.ldegree(y); j<=pol.degree(y); ++j) 00829 { 00830 for (int k=pol.ldegree(z); k<=pol.degree(z); ++k) 00831 { 00832 c = pol.coeff(x,i).coeff(y, j).coeff(z,k); 00833 if ( c != 0 ) cc.append(c); 00834 } 00835 } 00836 } 00837 return cc; 00838 } 00839 00840 00841 GiNaC::exvector coeff(GiNaC::ex pol) 00842 { 00843 using SyFi::x; 00844 using SyFi::y; 00845 using SyFi::z; 00846 00847 GiNaC::exvector cc; 00848 GiNaC::ex c, b; 00849 for (int i=pol.ldegree(x); i<=pol.degree(x); ++i) 00850 { 00851 for (int j=pol.ldegree(y); j<=pol.degree(y); ++j) 00852 { 00853 for (int k=pol.ldegree(z); k<=pol.degree(z); ++k) 00854 { 00855 c = pol.coeff(x,i).coeff(y, j).coeff(z,k); 00856 if ( c != 0 ) cc.insert(cc.begin(),c); 00857 } 00858 } 00859 } 00860 return cc; 00861 } 00862 00863 00864 GiNaC::exmap pol2basisandcoeff(GiNaC::ex e, GiNaC::ex s) 00865 { 00866 if (GiNaC::is_a<GiNaC::symbol>(s)) 00867 { 00868 GiNaC::symbol ss = GiNaC::ex_to<GiNaC::symbol>(s); 00869 e = expand(e); 00870 GiNaC::ex c; 00871 GiNaC::ex b; 00872 GiNaC::exmap map; 00873 for (int i=e.ldegree(ss); i<=e.degree(ss); ++i) 00874 { 00875 c = e.coeff(ss,i); 00876 b = pow(ss,i); 00877 map[b] = c; 00878 } 00879 return map; 00880 } 00881 else 00882 { 00883 throw(std::invalid_argument("The second argument must be a symbol.")); 00884 } 00885 } 00886 00887 00888 GiNaC::exmap pol2basisandcoeff(GiNaC::ex e) 00889 { 00890 using SyFi::x; 00891 using SyFi::y; 00892 using SyFi::z; 00893 00894 e = expand(e); 00895 GiNaC::ex c; 00896 GiNaC::ex b; 00897 GiNaC::exmap map; 00898 for (int i=e.ldegree(x); i<=e.degree(x); ++i) 00899 { 00900 for (int j=e.ldegree(y); j<=e.degree(y); ++j) 00901 { 00902 for (int k=e.ldegree(z); k<=e.degree(z); ++k) 00903 { 00904 c = e.coeff(x,i).coeff(y, j).coeff(z,k); 00905 b = pow(x,i)*pow(y,j)*pow(z,k); 00906 map[b] = c; 00907 } 00908 } 00909 } 00910 return map; 00911 } 00912 00913 00914 GiNaC::ex legendre1D(const GiNaC::symbol x, unsigned int n) 00915 { 00916 GiNaC::ex P; 00917 // Rodrigue's formula for Legendre polynomial of 1D 00918 // The interval [-1, 1] 00919 P=1/(pow(2,n)*GiNaC::factorial(n))*GiNaC::diff(GiNaC::pow((x*x-1),n),x,n); 00920 // ----------------- 00921 // The interval [0,1] 00922 // GiNaC::ex xx = 2*x - 1; 00923 // P=1/(pow(2,2*n)*GiNaC::factorial(n))*GiNaC::diff(GiNaC::pow((xx*xx-1),n),x,n); 00924 return P; 00925 } 00926 00927 00928 GiNaC::ex legendre(unsigned int order, unsigned int nsd, const string s) 00929 { 00930 using SyFi::x; 00931 using SyFi::y; 00932 using SyFi::z; 00933 00934 // The Legendre polynomials to be used in FiniteElement 00935 GiNaC::ex leg; 00936 GiNaC::ex A; 00937 GiNaC::lst basis; 00938 int dof; 00939 00940 GiNaC::ex b; 00941 00942 // 1D 00943 if(nsd == 1) 00944 { 00945 dof = order+1; 00946 A = get_symbolic_matrix(1,dof,s); 00947 int o=0; 00948 for(GiNaC::const_iterator i = A.begin(); i!=A.end(); ++i) 00949 { 00950 b= legendre1D(x,o); 00951 leg+= (*i)*b; 00952 basis.append(b); 00953 o++; 00954 } 00955 } 00956 // 2D 00957 /* 00958 else if(nsd == 2){ // NB: Only for tensor products on TRIANGLES (not boxes) 00959 / * 2D: structure of coefficients (a_i) 00960 * [ a_0 a_1 P_1(x) a_3 P_2(x) a_6 P_3(x) 00961 * [ a_2 P_1(y) a_4 P_1(x)*P_1(y) a_7 P_2(x)*P_1(y) 00962 * [ a_5 P_2(y) a_8 P_1(x)*P_2(y) 00963 * [ a_9 P_3(y) 00964 * / 00965 dof = (order+1)*(order+2)/2; 00966 A = get_symbolic_matrix(1,dof,s); 00967 size_t i=0; 00968 for (int o = 0; o <= order; o++) { 00969 for (int d = 0; d <= o; d++) { 00970 b = legendre1D(y,d)*legendre1D(x,o-d); 00971 leg += A.op(i)*b; 00972 basis.append(b); 00973 i++; 00974 00975 } 00976 } 00977 } 00978 */ 00979 else if(nsd == 2) // NB: Only for tensor products on rectangles 00980 { 00981 dof = (order+1)*(order+1); 00982 A = get_symbolic_matrix(1,dof,s); 00983 size_t i=0; 00984 for (unsigned int o = 0; o <= order; o++) 00985 { 00986 for (unsigned int d = 0; d <= order; d++) 00987 { 00988 b = legendre1D(y,d)*legendre1D(x,o); 00989 leg += A.op(i)*b; 00990 basis.append(b); 00991 i++; 00992 00993 } 00994 } 00995 } 00996 00997 /* tetrahedron 00998 else if(nsd==3){ 00999 dof = 0; 01000 for (int j=0; j<= order; j++) { 01001 dof += (j+1)*(j+2)/2; 01002 } 01003 A = get_symbolic_matrix(1, dof , s); 01004 01005 size_t i=0; 01006 for (int o = 0; o <= order; o++) { 01007 for (int d = 0; d <= o; d++) { 01008 for (int f = 0; f <= o; f++) { 01009 if ( o-d-f >= 0) { 01010 b = legendre1D(y,f)*legendre1D(z,d)*legendre1D(x,o-d-f); 01011 leg += A.op(i)*b; 01012 basis.append(b); 01013 i++; 01014 } 01015 } 01016 } 01017 } 01018 } 01019 */ 01020 01021 else if(nsd==3) 01022 { 01023 dof = (order+1)*(order+1)*(order+1); 01024 A = get_symbolic_matrix(1, dof , s); 01025 01026 size_t i=0; 01027 for (unsigned int o = 0; o <= order; o++) 01028 { 01029 for (unsigned int d = 0; d <= order; d++) 01030 { 01031 for (unsigned int f = 0; f <= order; f++) 01032 { 01033 b = legendre1D(y,f)*legendre1D(z,d)*legendre1D(x,o); 01034 leg += A.op(i)*b; 01035 basis.append(b); 01036 i++; 01037 } 01038 } 01039 } 01040 } 01041 return GiNaC::lst(leg,matrix_to_lst2(A), basis); 01042 } 01043 01044 01045 GiNaC::lst legendrev(unsigned int no_fields, unsigned int order, unsigned int nsd, const string a) 01046 { 01047 GiNaC::lst ret1; // contains the polynom 01048 GiNaC::lst ret2; // contains the coefficients 01049 GiNaC::lst ret3; // constains the basis functions 01050 GiNaC::lst basis_tmp; 01051 for (unsigned int i=1; i<= no_fields; i++) 01052 { 01053 GiNaC::lst basis; 01054 std::ostringstream s; 01055 s <<a<<""<<i<<"_"; 01056 GiNaC::ex polspace = legendre(order, nsd, s.str()); 01057 ret1.append(polspace.op(0)); 01058 ret2.append(polspace.op(1)); 01059 basis_tmp = GiNaC::ex_to<GiNaC::lst>(polspace.op(2)); 01060 for (GiNaC::lst::const_iterator basis_iterator = basis_tmp.begin(); 01061 basis_iterator != basis_tmp.end(); ++basis_iterator) 01062 { 01063 GiNaC::lst tmp_lst; 01064 for (unsigned int d=1; d<=no_fields; d++) tmp_lst.append(0); 01065 tmp_lst.let_op(i-1) = (*basis_iterator); 01066 ret3.append(tmp_lst); 01067 } 01068 } 01069 return GiNaC::lst(ret1,ret2,ret3); 01070 } 01071 01072 01073 bool compare(const ex & e, const string & s) 01074 { 01075 ostringstream ss; 01076 ss << e; 01077 return ss.str() == s; 01078 } 01079 01080 01081 void EQUAL_OR_DIE(const ex & e, const string & s) 01082 { 01083 if (!compare(e, s)) 01084 { 01085 ostringstream os; 01086 os << "ERROR: expression e: " <<e<<" is not equal to "<<s<<endl; 01087 throw runtime_error(os.str()); 01088 } 01089 } 01090 01091 01092 // =========== expression inspection utilities 01093 01094 class SymbolMapBuilderVisitor: 01095 public visitor, 01096 public symbol::visitor 01097 { 01098 public: 01099 SymbolMapBuilderVisitor(GiNaC::exmap & em): symbolmap(em) {} 01100 01101 GiNaC::exmap & symbolmap; 01102 01103 void visit(const symbol & s) 01104 { 01105 GiNaC::exmap::iterator it = symbolmap.find(s); 01106 if(it != symbolmap.end()) 01107 { 01108 it->second++; 01109 } 01110 else 01111 { 01112 //GiNaC::exmap::value_type p(ex(s), ex(s)); 01113 //symbolmap.insert(p); 01114 symbolmap[ex(s)] = ex(s); 01115 } 01116 } 01117 }; 01118 01119 class SymbolCounterVisitor: 01120 public visitor, 01121 public symbol::visitor, 01122 public basic::visitor 01123 { 01124 public: 01125 exhashmap<int> symbolcount; 01126 01127 void visit(const basic & s) 01128 { 01129 std::cout << "visiting basic " << std::endl; 01130 } 01131 01132 void visit(const symbol & s) 01133 { 01134 ex e = s; 01135 std::cout << "visiting symbol " << e << std::endl; 01136 exhashmap<int>::iterator it = symbolcount.find(s); 01137 if(it != symbolcount.end()) 01138 { 01139 std::cout << "found symbol " << e << std::endl; 01140 it->second++; 01141 } 01142 else 01143 { 01144 std::cout << "adding symbol " << e << std::endl; 01145 pair<ex,int> p; 01146 p.first = ex(s); 01147 p.second = 1; 01148 symbolcount.insert(p); 01149 } 01150 } 01151 }; 01152 01153 exhashmap<int> count_symbols(const ex & e) 01154 { 01155 SymbolCounterVisitor v; 01156 e.traverse(v); 01157 return v.symbolcount; 01158 } 01159 01160 01161 /* 01162 GiNaC::exvector get_symbols3(const GiNaC::ex & e) 01163 { 01164 // Implemented directly to avoid copying map: 01165 SymbolCounterVisitor v; 01166 e.traverse(v); 01167 exhashmap<int> & sc = v.symbolcount; 01168 01169 int i = 0; 01170 GiNaC::exvector l(sc.size()); 01171 for(exhashmap<int>::iterator it=sc.begin(); it!=sc.end(); it++) 01172 { 01173 //l.push_back(it->first); 01174 l[i] = it->first; 01175 i++; 01176 } 01177 return l; 01178 } 01179 01180 std::list<GiNaC::ex> get_symbols2(const ex & e) 01181 { 01182 // Implemented directly to avoid copying map: 01183 SymbolCounterVisitor v; 01184 e.traverse(v); 01185 exhashmap<int> & sc = v.symbolcount; 01186 01187 list<ex> l; 01188 for(exhashmap<int>::iterator it=sc.begin(); it!=sc.end(); it++) 01189 { 01190 l.push_back(it->first); 01191 } 01192 return l; 01193 } 01194 01195 void get_symbols(const ex & e, GiNaC::exmap & em) 01196 { 01197 SymbolMapBuilderVisitor v(em); 01198 e.traverse(v); 01199 } 01200 */ 01201 ex extract_symbols(const ex & e) 01202 { 01203 // Implemented directly to avoid copying map: 01204 SymbolCounterVisitor v; 01205 e.traverse(v); 01206 exhashmap<int> & sc = v.symbolcount; 01207 01208 lst l; 01209 for(exhashmap<int>::iterator it=sc.begin(); it!=sc.end(); it++) 01210 { 01211 l.append(it->first); 01212 std::cout << (it->first) << std::endl; 01213 } 01214 ex ret = l; 01215 return ret; 01216 } 01217 01218 01219 // Collect all symbols of an expression 01220 void collect_symbols(const GiNaC::ex & e, exset & v) 01221 { 01222 if (GiNaC::is_a<GiNaC::symbol>(e)) 01223 { 01224 v.insert(e); 01225 } 01226 else 01227 { 01228 for (size_t i=0; i<e.nops(); i++) 01229 { 01230 collect_symbols(e.op(i), v); 01231 } 01232 } 01233 } 01234 01235 01236 GiNaC::exvector collect_symbols(const GiNaC::ex & e) 01237 { 01238 exset s; 01239 collect_symbols(e, s); 01240 GiNaC::exvector v(s.size()); 01241 for(exset::iterator i=s.begin(); i!= s.end(); i++) 01242 { 01243 v.push_back(*i); 01244 } 01245 return v; 01246 } 01247 01248 01249 bool compare_archives(const string & first, const string & second, std::ostream & os) 01250 { 01251 bool ret = true; 01252 01253 // read both archives 01254 archive a1, a2; 01255 ifstream if1(first.c_str()), if2(second.c_str()); 01256 if1 >> a1; 01257 if2 >> a2; 01258 01259 // compare size 01260 int n = a1.num_expressions(); 01261 int n2 = a2.num_expressions(); 01262 if(n != n2) 01263 { 01264 os << "Archives " << first << " and " << second 01265 << " has a different number of expressions, " << n << " and " << n2 << "." << endl; 01266 os << "Comparing common expressions." << endl; 01267 ret = false; 01268 } 01269 01270 // iterate over all expressions in first archive 01271 ex e1,e2; 01272 for(int i=0; i<n; i++) 01273 { 01274 lst syms; 01275 string exname; 01276 01277 e1 = a1.unarchive_ex(syms, exname, i); 01278 01279 syms = ex_to<lst>(extract_symbols(e1)); 01280 // os << "Comparing " << exname << " with symbols " << syms << endl; 01281 01282 // is this in the second archive? 01283 try 01284 { 01285 e2 = a2.unarchive_ex(syms, exname.c_str()); 01286 01287 // got it, now compare 01288 bool isequal = is_zero(e1-e2); 01289 if(!isequal) 01290 { 01291 if(ret) 01292 { 01293 os << "Archives " << first << " and " << second 01294 << " are not equal, details follow:" << endl; 01295 } 01296 os << "Expression with name " << exname << " is not equal:" << endl; 01297 os << "First: " << endl << e1 << endl; 01298 os << "Second: " << endl << e2 << endl; 01299 ret = false; 01300 } 01301 } 01302 catch(...) 01303 { 01304 os << "Expression " << exname << " is missing from " << second << "." << endl; 01305 ret = false; 01306 } 01307 } 01308 01309 return ret; 01310 } 01311 01312 01313 class ExStatsVisitor: 01314 public visitor, 01315 public power::visitor, 01316 public mul::visitor, 01317 public add::visitor, 01318 public function::visitor 01319 { 01320 public: 01321 ExStats es; 01322 01323 void visit(const power & e) 01324 { 01325 //cout << "pow" << endl; 01326 ex b = e.op(1); 01327 ex c = b.evalf(); 01328 if( is_a<numeric>(c) ) 01329 { 01330 numeric nu = ex_to<numeric>(c); 01331 int n = (int) to_double(nu); 01332 if( is_zero(nu-n) ) 01333 { 01334 // is an integer power, interpret as a sequence of multiplications 01335 es.muls += n-1; 01336 es.flops += n-1; 01337 } 01338 } 01339 01340 es.pows++; 01341 } 01342 01343 void visit(const mul & e) 01344 { 01345 //cout << "mul" << endl; 01346 es.muls += e.nops()-1; 01347 es.flops += e.nops()-1; 01348 } 01349 01350 void visit(const add & e) 01351 { 01352 //cout << "add" << endl; 01353 es.adds += e.nops()-1; 01354 es.flops += e.nops()-1; 01355 } 01356 01357 void visit(const function & e) 01358 { 01359 //cout << "func" << endl; 01360 es.functions++; 01361 } 01362 }; 01363 01364 ExStats count_ops(const ex & e) 01365 { 01366 //cout << "count_ops " << e << endl; 01367 //cout << "is an add: " << GiNaC::is<GiNaC::add>(e) << endl; 01368 //cout << "is a mul: " << GiNaC::is<GiNaC::mul>(e) << endl; 01369 ExStatsVisitor v; 01370 e.traverse(v); 01371 return v.es; 01372 } 01373 01374 01375 // =========== expression manipulation utilities 01376 01377 // Returns in sel a list with symbol/expression pairs for all 01378 // power replacement variables introduced. Works best on 01379 // expressions in expanded form. 01380 ex replace_powers(const ex & ein, const list<symbol> & symbols, list<symexpair> & sel, const string & tmpsymbolprefix) 01381 { 01382 ex e = ein; 01383 // build power expressions 01384 list<symbol>::const_iterator it = symbols.begin(); 01385 for(; it != symbols.end(); it++) 01386 { 01387 int deg = e.degree(*it); 01388 if(deg > 0) 01389 { 01390 symbol sym = ex_to<symbol>(*it); 01391 string sname = tmpsymbolprefix + sym.get_name(); 01392 01393 // make list of new symbols 01394 vector<symbol> symbols(deg); 01395 symbols[0] = sym; 01396 for(int i=1; i<deg; i++) 01397 { 01398 symbols[i] = get_symbol( sname + int2string(i+1) ); 01399 sel.push_back(make_pair(symbols[i], symbols[i-1]*sym)); 01400 } 01401 01402 // with highest order first, subs in e 01403 ex prod = sym; 01404 for(int i=deg-1; i>=1; i--) 01405 { 01406 e = e.subs(power(sym,i+1) == symbols[i], subs_options::algebraic); 01407 } 01408 } 01409 } 01410 return e; 01411 } 01412 01413 01414 }; // namespace SyFi