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 <sstream> 00005 00006 #include "Polygon.h" 00007 #include "tools.h" 00008 #include "symbol_factory.h" 00009 00010 //using namespace std; 00011 00012 using std::string; 00013 using GiNaC::ex; 00014 using GiNaC::lst; 00015 using GiNaC::exvector; 00016 using GiNaC::ex_to; 00017 using GiNaC::is_a; 00018 using GiNaC::numeric; 00019 00020 namespace SyFi 00021 { 00022 00023 //--- Polygon ------------------------------------ 00024 // 00025 00026 Polygon::Polygon(const string & subscript_, unsigned int no_vert): 00027 subscript(subscript_), 00028 p(no_vert) 00029 { 00030 } 00031 00032 Polygon::Polygon(const Polygon& polygon) : 00033 subscript(polygon.str()), 00034 p(polygon.no_vertices()) 00035 { 00036 for (unsigned int i=0; i<p.size(); i++) 00037 { 00038 p[i] = polygon.vertex(i); 00039 } 00040 } 00041 00042 unsigned int Polygon::no_vertices() const 00043 { 00044 return p.size(); 00045 } 00046 00047 ex Polygon::vertex (unsigned int i) const 00048 { 00049 if ( i < 0 || i > p.size()-1 ) 00050 { 00051 throw std::out_of_range("Vertex index is out of range!"); 00052 } 00053 return p[i]; 00054 } 00055 00056 Line Polygon::line(unsigned int i) const 00057 { 00058 throw std::logic_error("line not implemented for this polygon subclass"); 00059 } 00060 00061 Triangle Polygon::triangle(unsigned int i) const 00062 { 00063 throw std::logic_error("triangle not implemented for this polygon subclass"); 00064 } 00065 00066 Rectangle Polygon::rectangle(unsigned int i) const 00067 { 00068 throw std::logic_error("rectangle not implemented for this polygon subclass"); 00069 } 00070 00071 //--- Line ------------------------------------ 00072 // 00073 00074 Line::Line(ex x0, ex x1, const string & subscript_): 00075 Polygon(subscript_, 2) 00076 { 00077 /* Lines on the form 00078 * x = x_0 + a_1 t 00079 * y = y_0 + a_2 t 00080 * z = z_0 + a_3 t 00081 */ 00082 00083 p[0] = x0; 00084 p[1] = x1; 00085 00086 //update internal structure 00087 if ( nsd == 1 ) 00088 { 00089 a_ = x1-x0; 00090 b_ = x0; 00091 } 00092 else if ( (GiNaC::is_a<lst>(x0) && GiNaC::is_a<lst>(x1)) && (x0.nops() == x1.nops()) ) 00093 { 00094 // TODO: support matrix in addition to lst? 00095 lst aa; 00096 lst bb; 00097 00098 // compute aa and bb 00099 for (unsigned int i=0; i<= x0.nops()-1; i++) 00100 { 00101 bb.append(x0.op(i)); 00102 aa.append(x1.op(i) - x0.op(i)); 00103 } 00104 00105 a_ = aa; 00106 b_ = bb; 00107 00108 } 00109 else 00110 { 00111 throw(std::logic_error("The points x0 and x1 needs to be of type lst if nsd > 1.")); 00112 } 00113 } 00114 00115 Line::Line(const Line& line): Polygon(line), a_(line.a_), b_(line.b_) 00116 { 00117 } 00118 00119 unsigned int Line::no_space_dim() const { return 1; } 00120 00121 ex Line::a() const 00122 { 00123 return a_; 00124 } 00125 00126 ex Line::b() const 00127 { 00128 return b_; 00129 } 00130 00131 ex Line::repr(Repr_format format) const 00132 { 00133 return repr(GiNaC::symbol("t"), format); 00134 } 00135 00136 ex Line::repr(ex t, Repr_format format) const 00137 { 00138 /* NOTE: use the same symbols in this 00139 * description as in the actual symbols 00140 * below and in the documentation. 00141 * 00142 * Lines on the form 00143 * x = x_0 + a_1 t, 00144 * y = y_0 + a_2 t, 00145 * z = z_0 + a_3 t. 00146 * for t in [0,1] 00147 */ 00148 00149 lst r; 00150 00151 if ( format == SUBS_PERFORMED ) 00152 { 00153 if (!GiNaC::is_a<lst>(a_)) 00154 { 00155 r = lst( x == b_ + a_*t); 00156 } 00157 else 00158 { 00159 if ( a_.nops() == 1) 00160 { 00161 r = lst( x == b_.op(0) + a_.op(0)*t); 00162 } // 2D 00163 else if ( a_.nops() == 2) 00164 { 00165 r = lst( x == b_.op(0) + a_.op(0)*t, 00166 y == b_.op(1) + a_.op(1)*t); 00167 // 3D 00168 } 00169 else if ( a_.nops() == 3) 00170 { 00171 r = lst( x == b_.op(0) + a_.op(0)*t, 00172 y == b_.op(1) + a_.op(1)*t, 00173 z == b_.op(2) + a_.op(2)*t); 00174 } 00175 } 00176 r.append(lst(t,0,1)); 00177 } 00178 else if ( format == SUBS_NOT_PERFORMED ) 00179 { 00180 00181 // NOTE: decide how the constant should be !! 00182 GiNaC::symbol a1("A"+subscript); 00183 GiNaC::symbol a2("C"+subscript); 00184 GiNaC::symbol a3("E"+subscript); 00185 GiNaC::symbol b1("B"+subscript); 00186 GiNaC::symbol b2("D"+subscript); 00187 GiNaC::symbol b3("F"+subscript); 00188 00189 // 2D 00190 if ( a_.nops() == 2) 00191 { 00192 r = lst( x == b1 + a1*t, 00193 y == b2 + a2*t); 00194 00195 r.append( a1 == a_.op(0)); 00196 r.append( b1 == b_.op(0)); 00197 r.append( a2 == a_.op(1)); 00198 r.append( b2 == b_.op(1)); 00199 // 3D 00200 } 00201 else if ( a_.nops() == 3) 00202 { 00203 r = lst( x == b1 + a1*t, 00204 y == b2 + a2*t, 00205 z == b3 + a3*t); 00206 00207 r.append( a1 == a_.op(0)); 00208 r.append( b1 == b_.op(0)); 00209 r.append( a2 == a_.op(1)); 00210 r.append( b2 == b_.op(1)); 00211 r.append( a3 == a_.op(2)); 00212 r.append( b3 == b_.op(2)); 00213 } 00214 r.append(lst(t,0,1)); 00215 } 00216 return r; 00217 } 00218 00219 const string Line::str() const 00220 { 00221 std::ostringstream s; 00222 // s <<"Line("<<p[0]<<","<<p[1]<<")"; 00223 // FIXME: would like to use the above code, but this causes a strange crash in Python 00224 s <<"Line"; 00225 return s.str(); 00226 } 00227 00228 ex Line::integrate(ex f, Repr_format format) 00229 { 00230 ex ret; 00231 00232 if ( format == SUBS_PERFORMED) 00233 { 00234 GiNaC::symbol t("t"); 00235 ex t_repr = repr(t); 00236 lst sub; 00237 int counter=0; 00238 // perform substitution 00239 if ( p[0].nops() == 3) 00240 { 00241 sub = lst(t_repr.op(0), t_repr.op(1), t_repr.op(2)); 00242 counter = 3; 00243 } 00244 else if ( p[0].nops() == 2) 00245 { 00246 sub = lst(t_repr.op(0), t_repr.op(1)); 00247 counter = 2; 00248 } 00249 else if ( p[0].nops() == 1) 00250 { 00251 sub = lst(t_repr.op(0)); 00252 counter = 1; 00253 } 00254 else if ( p[0].nops() == 0) 00255 { 00256 sub = lst(t_repr.op(0)); 00257 counter = 1; 00258 } 00259 00260 // compute D 00261 ex D; 00262 if ( p[0].nops() == 3) 00263 { 00264 D = pow(t_repr.op(0).rhs().coeff(t), 2) + 00265 pow(t_repr.op(1).rhs().coeff(t), 2) + 00266 pow(t_repr.op(2).rhs().coeff(t), 2); 00267 } 00268 else if ( p[0].nops() == 2) 00269 { 00270 D = pow(t_repr.op(0).rhs().coeff(t), 2) + 00271 pow(t_repr.op(1).rhs().coeff(t), 2); 00272 } 00273 else if ( p[0].nops() == 1) 00274 { 00275 D = pow(t_repr.op(0).rhs().coeff(t), 2); 00276 } 00277 else if ( p[0].nops() == 0) 00278 { 00279 D = pow(t_repr.op(0).rhs().coeff(t), 2); 00280 } 00281 00282 D = sqrt(D); 00283 00284 ex intf = f.subs(sub)*D; 00285 00286 intf = GiNaC::integral(t_repr.op(counter).op(0), t_repr.op(counter).op(1), t_repr.op(counter).op(2), intf); 00287 ret = eval_integ(intf); 00288 00289 } 00290 else if ( format == SUBS_NOT_PERFORMED ) 00291 { 00292 GiNaC::symbol t("t"); 00293 ex t_repr = repr(t); 00294 lst sub; 00295 GiNaC::symbol a("a"), b("b"), c("c"), D("D"); 00296 int counter=0; 00297 // perform substitution 00298 00299 if ( p[0].nops() == 3) 00300 { 00301 x == p[0].op(0) + a*t, 00302 // sub = lst(t_repr.op(0), t_repr.op(1), t_repr.op(2)); 00303 sub = lst( x == p[0].op(0) + a*t, 00304 y == p[0].op(1) + b*t, 00305 z == p[0].op(2) + c*t); 00306 counter = 3; 00307 } 00308 else if ( p[0].nops() == 2) 00309 { 00310 // sub = lst(t_repr.op(0), t_repr.op(1)); 00311 sub = lst( x == p[0].op(0) + a*t, 00312 y == p[0].op(1) + b*t); 00313 counter = 2; 00314 } 00315 00316 // compute D 00317 ex DD; 00318 if ( p[0].nops() == 3) 00319 { 00320 DD = pow(t_repr.op(0).rhs().coeff(t), 2) + 00321 pow(t_repr.op(1).rhs().coeff(t), 2) + 00322 pow(t_repr.op(2).rhs().coeff(t), 2); 00323 } 00324 else if ( p[0].nops() == 2) 00325 { 00326 DD = pow(t_repr.op(0).rhs().coeff(t), 2) + 00327 pow(t_repr.op(1).rhs().coeff(t), 2); 00328 } 00329 00330 DD = sqrt(DD); 00331 00332 ex intf = f.subs(sub); 00333 00334 intf = GiNaC::integral(t_repr.op(counter).op(0), t_repr.op(counter).op(1), t_repr.op(counter).op(2), intf); 00335 intf = eval_integ(intf); 00336 00337 ret = intf*D; 00338 00339 } 00340 else 00341 { 00342 throw std::runtime_error("Invalid format!"); 00343 } 00344 00345 return ret; 00346 } 00347 00348 Line* Line::copy() const 00349 { 00350 return new Line(*this); 00351 } 00352 00353 //--- ReferenceLine -------------------------------- 00354 // 00355 00356 ReferenceLine::ReferenceLine(const string & subscript_): 00357 Line(lst(0,0,0), lst(1,0,0), subscript_) 00358 { 00359 } 00360 00361 ReferenceLine::ReferenceLine(const ReferenceLine& line): Line(line) { } 00362 00363 ex ReferenceLine::repr(ex t, Repr_format format) const 00364 { 00365 lst r; 00366 r = lst( x == t, y == 0, z == 0, t, 0, 1); 00367 return r; 00368 } 00369 00370 const string ReferenceLine::str() const 00371 { 00372 std::ostringstream s; 00373 // s <<"ReferenceLine("<<p[0]<<","<<p[1]<<")"; 00374 s <<"ReferenceLine"; 00375 return s.str(); 00376 } 00377 00378 ex ReferenceLine::integrate(ex f, Repr_format formt) 00379 { 00380 ex intf = GiNaC::integral(x,0,1,f); 00381 intf = eval_integ(intf); 00382 return intf; 00383 } 00384 00385 ReferenceLine* ReferenceLine::copy() const 00386 { 00387 return new ReferenceLine(*this); 00388 } 00389 00390 //--- Triangle ------------------------------------ 00391 00392 Triangle::Triangle(const string & subscript_): 00393 Polygon(subscript_, 3) 00394 { 00395 // only used for the implementation of ReferenceTriangle 00396 } 00397 00398 Triangle::Triangle(ex x0, ex x1, ex x2, const string & subscript_): 00399 Polygon(subscript_, 3) 00400 { 00401 p[0] = x0; 00402 p[1] = x1; 00403 p[2] = x2; 00404 } 00405 00406 Triangle::Triangle(const Triangle& triangle): Polygon(triangle) { } 00407 00408 unsigned int Triangle::no_space_dim() const { return 2; } 00409 00410 Line Triangle::line(unsigned int i) const 00411 { 00412 00413 if ( i == 0) return Line(p[1],p[2], istr(subscript,i)); 00414 else if ( i == 1) return Line(p[0],p[2], istr(subscript,i)); 00415 else if ( i == 2) return Line(p[0],p[1], istr(subscript,i)); 00416 00417 throw std::out_of_range("Line index is out of range!"); 00418 } 00419 00420 ex Triangle::repr(Repr_format format) const 00421 { 00422 GiNaC::symbol r("r"), s("s"); 00423 GiNaC::symbol G("G"), H("H"); 00424 ex l1_repr = Line(vertex(0), vertex(1)).repr(r); 00425 ex l2_repr = Line(vertex(0), vertex(2)).repr(s); 00426 lst ret; 00427 //2D 00428 if ( p[0].nops() == 2) 00429 { 00430 ret = lst( x == l1_repr.op(0).rhs().coeff(r,0) 00431 + l1_repr.op(0).rhs().coeff(r,1)*r 00432 + l2_repr.op(0).rhs().coeff(s,1)*s, 00433 y == l1_repr.op(1).rhs().coeff(r,0) 00434 + l1_repr.op(1).rhs().coeff(r,1)*r 00435 + l2_repr.op(1).rhs().coeff(s,1)*s); 00436 00437 } 00438 else if ( p[0].nops() == 3) 00439 { 00440 00441 ret = lst( x == l1_repr.op(0).rhs().coeff(r,0) 00442 + l1_repr.op(0).rhs().coeff(r,1)*r 00443 + l2_repr.op(0).rhs().coeff(s,1)*s, 00444 y == l1_repr.op(1).rhs().coeff(r,0) 00445 + l1_repr.op(1).rhs().coeff(r,1)*r 00446 + l2_repr.op(1).rhs().coeff(s,1)*s, 00447 z == l1_repr.op(2).rhs().coeff(r,0) 00448 + l1_repr.op(2).rhs().coeff(r,1)*r 00449 + l2_repr.op(2).rhs().coeff(s,1)*s); 00450 00451 } 00452 00453 ret.append(lst(r, 0, 1)); 00454 ret.append(lst(s, 0, 1 - r)); 00455 00456 return ret; 00457 00458 } 00459 00460 const string Triangle::str() const 00461 { 00462 std::ostringstream s; 00463 // s <<"Triangle("<<p[0]<<","<<p[1]<<","<<p[2]<<")"; 00464 s <<"Triangle"; 00465 return s.str(); 00466 } 00467 00468 ex Triangle::integrate(ex func, Repr_format format) 00469 { 00470 ex ret; 00471 00472 if ( format == SUBS_PERFORMED ) 00473 { 00474 ex t_repr = repr(); 00475 lst sub; 00476 int counter=0; 00477 00478 // perform substitution 00479 if ( p[0].nops() == 3) 00480 { 00481 sub = lst(t_repr.op(0), t_repr.op(1), t_repr.op(2)); 00482 counter = 3; 00483 } 00484 else if ( p[0].nops() == 2) 00485 { 00486 sub = lst(t_repr.op(0), t_repr.op(1)); 00487 counter = 2; 00488 } 00489 ex intf = func.subs(sub); 00490 00491 // compute D 00492 ex D; 00493 if ( p[0].nops() == 3) 00494 { 00495 ex r = t_repr.op(3).op(0); 00496 ex s = t_repr.op(4).op(0); 00497 ex a = t_repr.op(0).rhs().coeff(r,1); 00498 ex b = t_repr.op(0).rhs().coeff(s,1); 00499 ex c = t_repr.op(1).rhs().coeff(r,1); 00500 ex d = t_repr.op(1).rhs().coeff(s,1); 00501 ex e = t_repr.op(2).rhs().coeff(r,1); 00502 ex f = t_repr.op(2).rhs().coeff(s,1); 00503 D = pow(c*f-d*e,2) + pow(a*f - b*e,2) + pow(a*d - b*c,2); 00504 D = sqrt(D); 00505 } 00506 else if ( p[0].nops() == 2) 00507 { 00508 ex r = t_repr.op(2).op(0); 00509 ex s = t_repr.op(3).op(0); 00510 ex a = t_repr.op(0).rhs().coeff(r,1); 00511 ex b = t_repr.op(0).rhs().coeff(s,1); 00512 ex c = t_repr.op(1).rhs().coeff(r,1); 00513 ex d = t_repr.op(1).rhs().coeff(s,1); 00514 D = abs(a*d-b*c); 00515 } 00516 00517 intf = intf*D; 00518 00519 counter++; 00520 intf = GiNaC::integral(t_repr.op(counter).op(0), t_repr.op(counter).op(1), t_repr.op(counter).op(2), intf); 00521 intf = eval_integ(intf); 00522 00523 counter--; 00524 intf = GiNaC::integral(t_repr.op(counter).op(0), t_repr.op(counter).op(1), t_repr.op(counter).op(2), intf); 00525 ret = eval_integ(intf); 00526 00527 } 00528 else if ( format == SUBS_NOT_PERFORMED ) 00529 { 00530 ex t_repr = repr(); 00531 lst sub; 00532 int counter=0; 00533 00534 GiNaC::symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f"), r("r"), s("s"), D("D"); 00535 00536 // perform substitution 00537 if ( p[0].nops() == 3) 00538 { 00539 sub = lst(x == p[0].op(0) + a*r + b*s, 00540 y == p[0].op(1) + c*r + d*s, 00541 z == p[0].op(2) + e*r + f*s); 00542 counter = 3; 00543 } 00544 else if ( p[0].nops() == 2) 00545 { 00546 sub = lst(t_repr.op(0), t_repr.op(1)); 00547 sub = lst(x == p[0].op(0) + a*r + b*s, 00548 y == p[0].op(1) + c*r + d*s); 00549 counter = 2; 00550 } 00551 00552 ex intf = func.subs(sub); 00553 00554 counter++; 00555 intf = GiNaC::integral(t_repr.op(counter).op(0), t_repr.op(counter).op(1), t_repr.op(counter).op(2), intf); 00556 intf = eval_integ(intf); 00557 00558 counter--; 00559 intf = GiNaC::integral(t_repr.op(counter).op(0), t_repr.op(counter).op(1), t_repr.op(counter).op(2), intf); 00560 intf = eval_integ(intf); 00561 00562 ret = intf*D; 00563 00564 } 00565 else 00566 { 00567 throw std::runtime_error("Invalid format!"); 00568 } 00569 00570 return ret; 00571 } 00572 00573 Triangle* Triangle::copy() const 00574 { 00575 return new Triangle(*this); 00576 } 00577 00578 //--- ReferenceTriangle ---------------------------- 00579 00580 ReferenceTriangle::ReferenceTriangle(const string & subscript_): 00581 Triangle(subscript_) 00582 { 00583 ex x0; 00584 ex x1; 00585 ex x2; 00586 if ( nsd == 3) 00587 { 00588 x0 = lst(0,0,0); 00589 x1 = lst(1,0,0); 00590 x2 = lst(0,1,0); 00591 } 00592 else if ( nsd == 2 ) 00593 { 00594 x0 = lst(0,0); 00595 x1 = lst(1,0); 00596 x2 = lst(0,1); 00597 } 00598 p[0] = x0; 00599 p[1] = x1; 00600 p[2] = x2; 00601 } 00602 00603 ReferenceTriangle::ReferenceTriangle(const ReferenceTriangle& triangle): Triangle(triangle) { } 00604 00605 const string ReferenceTriangle::str() const 00606 { 00607 std::ostringstream s; 00608 // s <<"ReferenceTriangle("<<p[0]<<","<<p[1]<<","<<p[2]<<")"; 00609 s <<"ReferenceTriangle"; 00610 return s.str(); 00611 } 00612 00613 ex ReferenceTriangle::integrate(ex f, Repr_format format) 00614 { 00615 ex ret = GiNaC::eval_integ( 00616 GiNaC::integral(y,0,1, 00617 GiNaC::eval_integ( 00618 GiNaC::integral(x,0,1-y,f)))); 00619 return ret; 00620 } 00621 00622 // ---------- Rectangle -------------------------------------- 00623 00624 Rectangle::Rectangle(const string & subscript_): 00625 Polygon(subscript_,4) 00626 { 00627 // only used for the implementation of ReferenceRectangle 00628 } 00629 00630 Rectangle::Rectangle(ex p0, ex p1, const string & subscript_): 00631 Polygon(subscript_,4) 00632 { 00633 ex x0; 00634 ex x1; 00635 ex x2; 00636 ex x3; 00637 00638 if (p0.nops() == 2 && p1.nops() == 2) 00639 { 00640 // Follows the UFC numbering convention for a reference quadrilateral: 00641 x0 = lst(p0.op(0), p0.op(1)); 00642 x1 = lst(p1.op(0), p0.op(1)); 00643 x2 = lst(p1.op(0), p1.op(1)); 00644 x3 = lst(p0.op(0), p1.op(1)); 00645 00646 } 00647 else if (p0.nops() == 3 && p1.nops() == 3) 00648 { 00649 if ( p0.op(0) == p1.op(0) ) 00650 { 00651 00652 x0 = lst(p0.op(0), p0.op(1), p0.op(2)); 00653 x1 = lst(p0.op(0), p1.op(1), p0.op(2)); 00654 x2 = lst(p0.op(0), p1.op(1), p1.op(2)); 00655 x3 = lst(p0.op(0), p0.op(1), p1.op(2)); 00656 00657 } 00658 else if ( p0.op(1) == p1.op(1) ) 00659 { 00660 00661 x0 = lst(p0.op(0), p0.op(1), p0.op(2)); 00662 x1 = lst(p1.op(0), p0.op(1), p0.op(2)); 00663 x2 = lst(p1.op(0), p0.op(1), p1.op(2)); 00664 x3 = lst(p0.op(0), p0.op(1), p1.op(2)); 00665 00666 } 00667 else if ( p0.op(2) == p1.op(2) ) 00668 { 00669 00670 x0 = lst(p0.op(0), p0.op(1), p0.op(2)); 00671 x1 = lst(p1.op(0), p0.op(1), p0.op(2)); 00672 x2 = lst(p1.op(0), p1.op(1), p0.op(2)); 00673 x3 = lst(p0.op(0), p1.op(1), p0.op(2)); 00674 } 00675 } 00676 else 00677 { 00678 throw(std::invalid_argument("The points p0 and p1 must be of dimention either 2 or 3.")); 00679 } 00680 00681 p[0] = x0; 00682 p[1] = x1; 00683 p[2] = x2; 00684 p[3] = x3; 00685 } 00686 00687 Rectangle::Rectangle(ex p0, ex p1, ex p2, ex p3, const string & subscript_): 00688 Polygon(subscript_,4) 00689 { 00690 p[0] = p0; 00691 p[1] = p1; 00692 p[2] = p2; 00693 p[3] = p3; 00694 } 00695 00696 Rectangle::Rectangle(const Rectangle& rectangle): Polygon(rectangle) { } 00697 00698 Rectangle* Rectangle::copy() const 00699 { 00700 return new Rectangle(*this); 00701 } 00702 00703 unsigned int Rectangle::no_space_dim() const 00704 { 00705 return 2; 00706 } 00707 00708 Line Rectangle::line(unsigned int i) const 00709 { 00710 if ( i == 0) return Line(p[0],p[1], istr(subscript,i)); 00711 else if ( i == 1) return Line(p[1],p[2], istr(subscript,i)); 00712 else if ( i == 2) return Line(p[2],p[3], istr(subscript,i)); 00713 else if ( i == 3) return Line(p[3],p[0], istr(subscript,i)); 00714 00715 throw std::out_of_range("Line index is out of range!"); 00716 } 00717 00718 ex Rectangle::repr(Repr_format format) const 00719 { 00720 lst ret; 00721 GiNaC::symbol r("r"), s("s"), t("t"); 00722 if ( p[0].nops() == 2 ) 00723 { 00724 ret.append( x == p[0].op(0) + r*( p[2].op(0) - p[0].op(0))); 00725 ret.append( y == p[0].op(1) + s*( p[2].op(1) - p[0].op(1))); 00726 ret.append( lst(r,0,1) ); 00727 ret.append( lst(s,0,1) ); 00728 } 00729 else if ( p[0].nops() == 3 ) 00730 { 00731 ret.append( x == p[0].op(0) + r*( p[2].op(0) - p[0].op(0))); 00732 ret.append( y == p[0].op(1) + s*( p[2].op(1) - p[0].op(1))); 00733 ret.append( z == p[0].op(2) + t*( p[2].op(2) - p[0].op(2))); 00734 ret.append( lst(r,0,1) ); 00735 ret.append( lst(s,0,1) ); 00736 ret.append( lst(t,0,1) ); 00737 } 00738 return ret; 00739 } 00740 00741 const string Rectangle::str() const 00742 { 00743 std::ostringstream s; 00744 // s <<"Rectangle("<<p[0]<<","<<p[1]<<","<<p[2]<<","<<p[3]<<")"; 00745 s <<"Rectangle"; 00746 return s.str(); 00747 } 00748 00749 ex Rectangle::integrate(ex func, Repr_format format) 00750 { 00751 00752 unsigned int counter = 0; 00753 ex s_repr = repr(); 00754 00755 lst sub; 00756 // perform substitution 00757 if ( p[0].nops() == 3) 00758 { 00759 sub = lst(s_repr.op(0), s_repr.op(1), s_repr.op(2)); 00760 counter = 3; 00761 } 00762 else if ( p[0].nops() == 2) 00763 { 00764 sub = lst(s_repr.op(0), s_repr.op(1)); 00765 counter = 2; 00766 } 00767 00768 ex D; 00769 if ( p[0].nops() == 2) 00770 { 00771 D = ( p[2].op(0) - p[0].op(0)) 00772 *( p[2].op(1) - p[0].op(1)); 00773 } 00774 else if ( p[0].nops() == 3 ) 00775 { 00776 00777 if ( p[2].op(0) == p[0].op(0) ) 00778 { 00779 D = ( p[2].op(1) - p[0].op(1)) 00780 *( p[2].op(2) - p[0].op(2)); 00781 } 00782 else if ( p[2].op(1) == p[0].op(1) ) 00783 { 00784 D = ( p[2].op(0) - p[0].op(0)) 00785 *( p[2].op(2) - p[0].op(2)); 00786 } 00787 else if ( p[2].op(2) == p[0].op(2) ) 00788 { 00789 D = ( p[2].op(0) - p[0].op(0)) 00790 *( p[2].op(1) - p[0].op(1)); 00791 } 00792 } 00793 00794 ex intf = func.subs(sub); 00795 00796 intf = intf*D; 00797 00798 intf = GiNaC::integral(s_repr.op(counter).op(0), s_repr.op(counter).op(1), s_repr.op(counter).op(2), intf); 00799 intf = eval_integ(intf); 00800 00801 counter++; 00802 intf = GiNaC::integral(s_repr.op(counter).op(0), s_repr.op(counter).op(1), s_repr.op(counter).op(2), intf); 00803 intf = eval_integ(intf); 00804 00805 counter++; 00806 if ( counter < s_repr.nops() ) 00807 { 00808 intf = GiNaC::integral(s_repr.op(counter).op(0), s_repr.op(counter).op(1), s_repr.op(counter).op(2), intf); 00809 intf = eval_integ(intf); 00810 } 00811 00812 return intf; 00813 } 00814 00815 // ---------- ReferenceRectangle -------------------------------------- 00816 // 00817 00818 ReferenceRectangle::ReferenceRectangle(const string & subscript_): 00819 Rectangle(subscript_) 00820 { 00821 p[0] = lst(0,0); 00822 p[1] = lst(1,0); 00823 p[2] = lst(1,1); 00824 p[3] = lst(0,1); 00825 } 00826 00827 ReferenceRectangle::ReferenceRectangle(const ReferenceRectangle& rectangle): Rectangle(rectangle) { } 00828 00829 ReferenceRectangle* ReferenceRectangle::copy() const 00830 { 00831 return new ReferenceRectangle(*this); 00832 } 00833 00834 const string ReferenceRectangle::str() const 00835 { 00836 std::ostringstream s; 00837 // s <<"ReferenceRectangle("<<p[0]<<","<<p[1]<<","<<p[2]<<","<<p[3]<<")"; 00838 s <<"ReferenceRectangle"; 00839 return s.str(); 00840 } 00841 00842 ReferenceTriangle* ReferenceTriangle::copy() const 00843 { 00844 return new ReferenceTriangle(*this); 00845 } 00846 00847 //--- Tetrahedron ------------------------------------ 00848 00849 Tetrahedron::Tetrahedron(ex x0, ex x1, ex x2, ex x3, const string & subscript_): 00850 Polygon(subscript_,4) 00851 { 00852 p[0] = x0; 00853 p[1] = x1; 00854 p[2] = x2; 00855 p[3] = x3; 00856 } 00857 00858 Tetrahedron::Tetrahedron(const Tetrahedron& tetrahedron): Polygon(tetrahedron) { } 00859 00860 unsigned int Tetrahedron::no_space_dim() const 00861 { 00862 return 3; 00863 } 00864 00865 Line Tetrahedron::line(unsigned int i) const 00866 { 00867 int i0, i1; 00868 switch(i) 00869 { 00870 case 0: i0 = 0; i1 = 1; break; 00871 case 1: i0 = 0; i1 = 2; break; 00872 case 2: i0 = 0; i1 = 3; break; 00873 case 3: i0 = 1; i1 = 2; break; 00874 case 4: i0 = 1; i1 = 3; break; 00875 case 5: i0 = 2; i1 = 3; break; 00876 default: 00877 throw std::out_of_range("Line index is out of range!"); 00878 } 00879 return Line(p[i0], p[i1], istr(subscript,i)); 00880 } 00881 00882 Triangle Tetrahedron::triangle(unsigned int i) const 00883 { 00884 00885 if ( i == 0 ) 00886 { 00887 return Triangle(p[1], p[2], p[3], istr(subscript,i)); 00888 } 00889 else if ( i == 1) 00890 { 00891 return Triangle(p[0], p[2], p[3], istr(subscript,i)); 00892 } 00893 else if ( i == 2) 00894 { 00895 return Triangle(p[0], p[1], p[3], istr(subscript,i)); 00896 } 00897 else if ( i == 3) 00898 { 00899 return Triangle(p[0], p[1], p[2], istr(subscript,i)); 00900 } 00901 00902 throw std::out_of_range("Face index is out of range!"); 00903 00904 } 00905 00906 ex Tetrahedron::repr(Repr_format format) const 00907 { 00908 GiNaC::symbol r("r"), s("s"), t("t"); 00909 ex l1_repr = Line(vertex(0), vertex(1)).repr(r); 00910 ex l2_repr = Line(vertex(0), vertex(2)).repr(s); 00911 ex l3_repr = Line(vertex(0), vertex(3)).repr(t); 00912 lst ret; 00913 00914 ret = lst( 00915 x == l1_repr.op(0).rhs().coeff(r,0) + l1_repr.op(0).rhs().coeff(r,1)*r 00916 + l2_repr.op(0).rhs().coeff(s,1)*s + l3_repr.op(0).rhs().coeff(t,1)*t, 00917 y == l1_repr.op(1).rhs().coeff(r,0) + l1_repr.op(1).rhs().coeff(r,1)*r 00918 + l2_repr.op(1).rhs().coeff(s,1)*s + l3_repr.op(1).rhs().coeff(t,1)*t, 00919 z == l1_repr.op(2).rhs().coeff(r,0) + l1_repr.op(2).rhs().coeff(r,1)*r 00920 + l2_repr.op(2).rhs().coeff(s,1)*s + l3_repr.op(2).rhs().coeff(t,1)*t); 00921 00922 ret.append(lst(r, 0, 1)); 00923 ret.append(lst(s, 0, 1 - r)); 00924 ret.append(lst(t, 0, 1 - r - s)); 00925 00926 return ret; 00927 } 00928 00929 const string Tetrahedron::str() const 00930 { 00931 std::ostringstream s; 00932 // s <<"Tetrahedron("<<p[0]<<","<<p[1]<<","<<p[2]<<","<<p[3]<<")"; 00933 s <<"Tetrahedron"; 00934 return s.str(); 00935 } 00936 00937 ex Tetrahedron::integrate(ex func, Repr_format format) 00938 { 00939 ex ret; 00940 00941 if ( format == SUBS_PERFORMED ) 00942 { 00943 ex t_repr = repr(); 00944 // std::cout <<"t "<<t_repr<<std::endl; 00945 00946 //perform substitution 00947 lst sub = lst(t_repr.op(0), t_repr.op(1), t_repr.op(2)); 00948 ex intf = func.subs(sub); 00949 00950 // compute D 00951 ex D; 00952 ex r = t_repr.op(3).op(0); 00953 ex s = t_repr.op(4).op(0); 00954 ex t = t_repr.op(5).op(0); 00955 ex a = t_repr.op(0).rhs().coeff(r,1); 00956 ex b = t_repr.op(0).rhs().coeff(s,1); 00957 ex c = t_repr.op(0).rhs().coeff(t,1); 00958 ex d = t_repr.op(1).rhs().coeff(r,1); 00959 ex e = t_repr.op(1).rhs().coeff(s,1); 00960 ex f = t_repr.op(1).rhs().coeff(t,1); 00961 ex g = t_repr.op(2).rhs().coeff(r,1); 00962 ex h = t_repr.op(2).rhs().coeff(s,1); 00963 ex k = t_repr.op(2).rhs().coeff(t,1); 00964 00965 D = a*(e*k-f*h) - b*(d*k-f*g) + c*(d*h - g*e); 00966 00967 intf = intf*D; 00968 00969 intf = GiNaC::integral(t_repr.op(5).op(0), t_repr.op(5).op(1), t_repr.op(5).op(2), intf); 00970 intf = GiNaC::eval_integ(intf); 00971 00972 intf = GiNaC::integral(t_repr.op(4).op(0), t_repr.op(4).op(1), t_repr.op(4).op(2), intf); 00973 intf = GiNaC::eval_integ(intf); 00974 00975 intf = GiNaC::integral(t_repr.op(3).op(0), t_repr.op(3).op(1), t_repr.op(3).op(2), intf); 00976 ret = GiNaC::eval_integ(intf); 00977 00978 } 00979 else if ( format == SUBS_NOT_PERFORMED ) 00980 { 00981 ex t_repr = repr(); 00982 00983 GiNaC::symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f"), g("g"), h("h"), k("k"), D("D"); 00984 ex r = t_repr.op(3).op(0); 00985 ex s = t_repr.op(4).op(0); 00986 ex t = t_repr.op(5).op(0); 00987 00988 //perform substitution 00989 // lst sub = lst(t_repr.op(0), t_repr.op(1), t_repr.op(2)); 00990 ex sub = lst( 00991 x == p[0].op(0) + a*r + b*s + c*t, 00992 y == p[0].op(1) + d*r + e*s + f*t, 00993 z == p[0].op(2) + g*r + h*s + k*t); 00994 00995 ex intf = func.subs(sub); 00996 00997 intf = GiNaC::integral(t_repr.op(5).op(0), t_repr.op(5).op(1), t_repr.op(5).op(2), intf); 00998 intf = GiNaC::eval_integ(intf); 00999 01000 intf = GiNaC::integral(t_repr.op(4).op(0), t_repr.op(4).op(1), t_repr.op(4).op(2), intf); 01001 intf = GiNaC::eval_integ(intf); 01002 01003 intf = GiNaC::integral(t_repr.op(3).op(0), t_repr.op(3).op(1), t_repr.op(3).op(2), intf); 01004 intf = GiNaC::eval_integ(intf); 01005 01006 ret = intf*D; 01007 01008 } 01009 else 01010 { 01011 throw std::runtime_error("Invalid format."); 01012 } 01013 01014 return ret; 01015 } 01016 01017 Tetrahedron* Tetrahedron::copy() const 01018 { 01019 return new Tetrahedron(*this); 01020 } 01021 01022 //--- ReferenceTetrahedron --------------------------- 01023 // 01024 01025 ReferenceTetrahedron::ReferenceTetrahedron(const string & subscript_): 01026 Tetrahedron(lst(0, 0, 0), lst(1, 0, 0), lst(0, 1, 0), lst(0, 0, 1), subscript_) 01027 { 01028 } 01029 01030 ReferenceTetrahedron::ReferenceTetrahedron(const ReferenceTetrahedron& tetrahedron): Tetrahedron(tetrahedron) { } 01031 01032 const string ReferenceTetrahedron::str() const 01033 { 01034 std::ostringstream s; 01035 // s <<"ReferenceTetrahedron("<<p[0]<<","<<p[1]<<","<<p[2]<<","<<p[3]<<")"; 01036 s <<"ReferenceTetrahedron"; 01037 return s.str(); 01038 } 01039 01040 ex ReferenceTetrahedron::integrate(ex f, Repr_format format) 01041 { 01042 01043 ex intf = GiNaC::integral(x,0,1-y-z,f); 01044 intf = GiNaC::eval_integ(intf); 01045 01046 intf = GiNaC::integral(y,0,1-z, intf); 01047 intf = GiNaC::eval_integ(intf); 01048 01049 intf = GiNaC::integral(z,0,1, intf); 01050 intf = GiNaC::eval_integ(intf); 01051 01052 return intf; 01053 } 01054 01055 ReferenceTetrahedron* ReferenceTetrahedron::copy() const 01056 { 01057 return new ReferenceTetrahedron(*this); 01058 } 01059 01060 // ---------- Box -------------------------------------- 01061 // 01062 01063 Box::Box(ex p0, ex p1, const string & subscript_): 01064 Polygon(subscript_, 8) 01065 { 01066 // Numbered by the UFC convention: 01067 p[0] = lst(p0.op(0), p0.op(1), p0.op(2)); 01068 p[1] = lst(p1.op(0), p0.op(1), p0.op(2)); 01069 p[2] = lst(p1.op(0), p1.op(1), p0.op(2)); 01070 p[3] = lst(p0.op(0), p1.op(1), p0.op(2)); 01071 01072 p[4] = lst(p0.op(0), p0.op(1), p1.op(2)); 01073 p[5] = lst(p1.op(0), p0.op(1), p1.op(2)); 01074 p[6] = lst(p1.op(0), p1.op(1), p1.op(2)); 01075 p[7] = lst(p0.op(0), p1.op(1), p1.op(2)); 01076 } 01077 01078 Box::Box(ex p0, ex p1, ex p2, ex p3, ex p4, ex p5, ex p6, ex p7, const string & subscript_): 01079 Polygon(subscript_, 8) 01080 { 01081 p[0] = p0; 01082 p[1] = p1; 01083 p[2] = p2; 01084 p[3] = p3; 01085 01086 p[4] = p4; 01087 p[5] = p5; 01088 p[6] = p6; 01089 p[7] = p7; 01090 } 01091 01092 Box::Box(const Box& box): 01093 Polygon(box) 01094 { 01095 } 01096 01097 unsigned int Box::no_space_dim() const 01098 { 01099 return 3; 01100 } 01101 01102 // Numbered by the UFC convention: 01103 Line Box::line(unsigned int i) const 01104 { 01105 int i0, i1; 01106 switch(i) 01107 { 01108 case 0: i0 = 6; i1 = 7; break; 01109 case 1: i0 = 5; i1 = 6; break; 01110 case 2: i0 = 4; i1 = 7; break; 01111 case 3: i0 = 4; i1 = 5; break; 01112 case 4: i0 = 3; i1 = 7; break; 01113 case 5: i0 = 2; i1 = 6; break; 01114 case 6: i0 = 2; i1 = 3; break; 01115 case 7: i0 = 1; i1 = 5; break; 01116 case 8: i0 = 1; i1 = 2; break; 01117 case 9: i0 = 0; i1 = 4; break; 01118 case 10: i0 = 0; i1 = 3; break; 01119 case 11: i0 = 0; i1 = 1; break; 01120 default: 01121 throw std::out_of_range("Line index is out of range!"); 01122 } 01123 return Line(p[i0], p[i1], istr(subscript,i)); 01124 } 01125 01126 // Numbered by the UFC convention. 01127 Rectangle Box::rectangle(unsigned int i) const 01128 { 01129 switch(i) 01130 { 01131 case 0: return Rectangle(p[4], p[6], istr(subscript,i)); 01132 case 1: return Rectangle(p[2], p[7], istr(subscript,i)); 01133 case 2: return Rectangle(p[1], p[6], istr(subscript,i)); 01134 case 3: return Rectangle(p[0], p[7], istr(subscript,i)); 01135 case 4: return Rectangle(p[0], p[5], istr(subscript,i)); 01136 case 5: return Rectangle(p[0], p[2], istr(subscript,i)); 01137 } 01138 throw std::out_of_range("Rectangle index is out of range!"); 01139 } 01140 01141 ex Box::repr(Repr_format format) const 01142 { 01143 lst ret; 01144 GiNaC::symbol r("r"), s("s"), t("t"); 01145 ret.append( x == p[0].op(0) + r * (p[6].op(0) - p[0].op(0)) ); 01146 ret.append( y == p[0].op(1) + s * (p[6].op(1) - p[0].op(1)) ); 01147 ret.append( z == p[0].op(2) + t * (p[6].op(2) - p[0].op(2)) ); 01148 ret.append( lst(r,0,1) ); 01149 ret.append( lst(s,0,1) ); 01150 ret.append( lst(t,0,1) ); 01151 return ret; 01152 } 01153 01154 const string Box::str() const 01155 { 01156 std::ostringstream s; 01157 // s <<"Box("<<p[0]<<","<<p[1]<<","<<p[2]<<","<<p[3]<<","<<p[4]<<","<<p[5]<<","<<p[6]<<","<<p[7]<<")"; 01158 s <<"Box"; 01159 return s.str(); 01160 } 01161 01162 ex Box::integrate(ex func, Repr_format format) 01163 { 01164 01165 ex b_repr = repr(); 01166 01167 lst sub; 01168 sub = lst(b_repr.op(0), b_repr.op(1), b_repr.op(2)); 01169 ex intf = func.subs(sub); 01170 01171 ex D = ( p[6].op(0) - p[0].op(0) ) 01172 * ( p[6].op(1) - p[0].op(1) ) 01173 * ( p[6].op(2) - p[0].op(2) ); 01174 01175 intf = intf*D; 01176 01177 intf = GiNaC::integral(b_repr.op(3).op(0), b_repr.op(3).op(1), b_repr.op(3).op(2), intf); 01178 intf = eval_integ(intf); 01179 01180 intf = GiNaC::integral(b_repr.op(4).op(0), b_repr.op(4).op(1), b_repr.op(4).op(2), intf); 01181 intf = eval_integ(intf); 01182 01183 intf = GiNaC::integral(b_repr.op(5).op(0), b_repr.op(5).op(1), b_repr.op(5).op(2), intf); 01184 intf = eval_integ(intf); 01185 01186 return intf; 01187 } 01188 01189 Box* Box::copy() const 01190 { 01191 return new Box(*this); 01192 } 01193 01194 // ---------- ReferenceBox -------------------------------------- 01195 01196 ReferenceBox::ReferenceBox(const string & subscript_): 01197 Box(lst(0,0,0), lst(1,1,1), subscript_) 01198 { 01199 } 01200 01201 ReferenceBox::ReferenceBox(const ReferenceBox& box): 01202 Box(box) 01203 { 01204 } 01205 01206 const string ReferenceBox::str() const 01207 { 01208 std::ostringstream s; 01209 // s <<"ReferenceBox("<<p[0]<<","<<p[1]<<","<<p[2]<<","<<p[3]<<","<<p[4]<<","<<p[5]<<","<<p[6]<<","<<p[7]<<")"; 01210 s <<"ReferenceBox"; 01211 return s.str(); 01212 } 01213 01214 ReferenceBox* ReferenceBox::copy() const 01215 { 01216 return new ReferenceBox(*this); 01217 } 01218 01219 // ---------- General Simplex -------------------------------------- 01220 01221 Simplex::Simplex(GiNaC::lst vertices, const string & subscript_): 01222 Polygon(subscript_,vertices.nops()) 01223 { 01224 //FIXME does this work in 1D ? Need to check! 01225 unsigned int space_dim = vertices.op(0).nops(); 01226 for (unsigned int i=0; i<vertices.nops(); i++) 01227 { 01228 if (space_dim != vertices.op(i).nops()) 01229 { 01230 p.clear(); 01231 throw std::logic_error("Vertex dimension must be the same for all points!"); 01232 } 01233 p[i] = vertices.op(i); 01234 } 01235 } 01236 01237 Simplex::Simplex(const Simplex& simplex): 01238 Polygon(simplex) 01239 { 01240 } 01241 01242 unsigned int Simplex::no_space_dim() const 01243 { 01244 return p[0].nops()-1; 01245 } 01246 01247 ex Simplex::repr(Repr_format format) const 01248 { 01249 unsigned int nsd = vertex(0).nops(); 01250 unsigned int no_lines = no_vertices()-1; 01251 01252 ex r = get_symbolic_vector(nsd, "r"); 01253 ex x = get_symbolic_vector(nsd, "x"); 01254 ex ri; 01255 lst lines; 01256 for (unsigned int i=0; i< no_vertices()-1; i++) 01257 { 01258 ri = r.op(i); 01259 lst line_i_repr; 01260 for (unsigned int d=0; d< nsd; d++) 01261 { 01262 line_i_repr.append(x.op(d) == (vertex(i+1).op(d) - vertex(0).op(d))*ri + vertex(0).op(d)); 01263 } 01264 line_i_repr.append(lst(ri, 0, 1)); 01265 01266 lines.append(line_i_repr); 01267 } 01268 01269 lst ret; 01270 for (unsigned int i=0; i < nsd; i++) 01271 { 01272 ri = r.op(i); 01273 GiNaC::ex xi_expr; 01274 GiNaC::ex rhs = lines.op(0).op(i).rhs().coeff(ri,0); 01275 for (unsigned int l=0; l < no_lines; l++) 01276 { 01277 // xi_expr2 == xi_expr.lhs() == xi_expr.rhs() + lines.op(l).op(i).rhs().coeff(ri,1)*ri; 01278 rhs += lines.op(l).op(i).rhs().coeff(ri,1)*ri; 01279 01280 } 01281 xi_expr = x.op(i) == rhs; 01282 ret.append(xi_expr); 01283 } 01284 01285 GiNaC::ex limit=1; 01286 for (unsigned int i=0; i< no_lines; i++) 01287 { 01288 ri = r.op(i); 01289 ret.append(lst(ri, 0, limit)); 01290 limit -= ri; 01291 } 01292 01293 return ret; 01294 } 01295 01296 const string Simplex::str() const 01297 { 01298 std::ostringstream s; 01299 /* 01300 s <<"Simplex("; 01301 for (int i=0; i<p.size()-1; i++) { 01302 s << p[i]<<","; 01303 } 01304 s << p[p.size()-1]<<")"; 01305 */ 01306 s <<"Simplex"; 01307 return s.str(); 01308 } 01309 01310 ex Simplex::integrate(ex func, Repr_format format) 01311 { 01312 ex ret; 01313 01314 unsigned int nsd = vertex(0).nops(); 01315 01316 ex s_repr = repr(); 01317 lst subs; 01318 for (unsigned int i=0; i< nsd; i++) 01319 { 01320 subs.append(s_repr.op(i)); 01321 } 01322 01323 ex intf = func.subs(subs); 01324 01325 for (unsigned int i=s_repr.nops()-1; i>=nsd; i--) 01326 { 01327 intf = GiNaC::integral(s_repr.op(i).op(0), s_repr.op(i).op(1), s_repr.op(i).op(2), intf); 01328 intf = GiNaC::eval_integ(intf); 01329 } 01330 01331 return intf; 01332 } 01333 01334 Simplex Simplex::sub_simplex(unsigned int i) 01335 { 01336 if ( i < 0 || i >= p.size()) 01337 { 01338 throw std::out_of_range("Can only create subsimplices between 0 and the number of vertices-1."); 01339 } 01340 lst v; 01341 for (unsigned int k=0; k< p.size(); k++) 01342 { 01343 if ( k != i) 01344 { 01345 v.append(p[k]); 01346 } 01347 } 01348 return Simplex(v, istr(subscript, i)); 01349 } 01350 01351 Simplex* Simplex::copy() const 01352 { 01353 return new Simplex(*this); 01354 } 01355 01356 //-- tools for Polygons ------------------------------------------------------ 01357 01358 lst bezier_ordinates(Tetrahedron& tetrahedra, unsigned int d) 01359 { 01360 01361 //FIXME: ugly conversion to matrix 01362 01363 lst ret; 01364 ex V1 = tetrahedra.vertex(0); 01365 ex V2 = tetrahedra.vertex(1); 01366 ex V3 = tetrahedra.vertex(2); 01367 ex V4 = tetrahedra.vertex(3); 01368 01369 lst V1l = ex_to<lst>(V1); 01370 lst V2l = ex_to<lst>(V2); 01371 lst V3l = ex_to<lst>(V3); 01372 lst V4l = ex_to<lst>(V4); 01373 01374 ex V1m = lst_to_matrix2(V1l); 01375 ex V2m = lst_to_matrix2(V2l); 01376 ex V3m = lst_to_matrix2(V3l); 01377 ex V4m = lst_to_matrix2(V4l); 01378 01379 int l; 01380 for (unsigned int i=0; i<= d; i++) 01381 { 01382 for (unsigned int j=0; j<= d; j++) 01383 { 01384 for (unsigned int k=0; k<= d; k++) 01385 { 01386 if ( d - i - j -k >= 0 ) 01387 { 01388 l= d - i - j -k; 01389 ex sum = (l*V1m + k*V2m + j*V3m + i*V4m)/d; 01390 ret.append(matrix_to_lst2(sum.evalm())); 01391 } 01392 } 01393 } 01394 } 01395 // FIXME how should these be sorted ????? 01396 // ret = ret.sort(); 01397 return ret; 01398 } 01399 01400 lst interior_coordinates(Tetrahedron& tetrahedra, unsigned int d) 01401 { 01402 01403 //FIXME: ugly conversion to matrix 01404 d = d+4; 01405 01406 lst ret; 01407 ex V1 = tetrahedra.vertex(0); 01408 ex V2 = tetrahedra.vertex(1); 01409 ex V3 = tetrahedra.vertex(2); 01410 ex V4 = tetrahedra.vertex(3); 01411 01412 lst V1l = ex_to<lst>(V1); 01413 lst V2l = ex_to<lst>(V2); 01414 lst V3l = ex_to<lst>(V3); 01415 lst V4l = ex_to<lst>(V4); 01416 01417 ex V1m = lst_to_matrix2(V1l); 01418 ex V2m = lst_to_matrix2(V2l); 01419 ex V3m = lst_to_matrix2(V3l); 01420 ex V4m = lst_to_matrix2(V4l); 01421 01422 int l; 01423 for (unsigned int i=1; i< d; i++) 01424 { 01425 for (unsigned int j=1; j< d; j++) 01426 { 01427 for (unsigned int k=1; k< d; k++) 01428 { 01429 if ( d - i - j -k >= 1 ) 01430 { 01431 l= d - i - j -k; 01432 ex sum = (l*V1m + k*V2m + j*V3m + i*V4m)/d; 01433 ret.append(matrix_to_lst2(sum.evalm())); 01434 } 01435 } 01436 } 01437 } 01438 // FIXME how should these be sorted ????? 01439 // ret = ret.sort(); 01440 return ret; 01441 } 01442 01443 lst bezier_ordinates(Triangle& triangle, unsigned int d) 01444 { 01445 01446 //FIXME: ugly conversion to matrix 01447 01448 lst ret; 01449 ex V1 = triangle.vertex(0); 01450 ex V2 = triangle.vertex(1); 01451 ex V3 = triangle.vertex(2); 01452 01453 lst V1l = ex_to<lst>(V1); 01454 lst V2l = ex_to<lst>(V2); 01455 lst V3l = ex_to<lst>(V3); 01456 01457 ex V1m = lst_to_matrix2(V1l); 01458 ex V2m = lst_to_matrix2(V2l); 01459 ex V3m = lst_to_matrix2(V3l); 01460 01461 int k; 01462 for (unsigned int i=0; i <= d; i++) 01463 { 01464 for (unsigned int j=0; j <= d; j++) 01465 { 01466 if ( int(d) - int(i) - int(j) >= 0 ) 01467 { 01468 k = d - i - j; 01469 ex sum = (k*V1m + j*V2m + i*V3m)/d; 01470 ret.append(matrix_to_lst2(sum.evalm())); 01471 } 01472 } 01473 } 01474 // FIXME how should these be sorted ????? 01475 // ret = ret.sort(); 01476 return ret; 01477 } 01478 01479 lst interior_coordinates(Triangle& triangle, unsigned int d) 01480 { 01481 01482 //FIXME: ugly conversion to matrix 01483 // 01484 d=d+3; 01485 01486 lst ret; 01487 ex V1 = triangle.vertex(0); 01488 ex V2 = triangle.vertex(1); 01489 ex V3 = triangle.vertex(2); 01490 01491 lst V1l = ex_to<lst>(V1); 01492 lst V2l = ex_to<lst>(V2); 01493 lst V3l = ex_to<lst>(V3); 01494 01495 ex V1m = lst_to_matrix2(V1l); 01496 ex V2m = lst_to_matrix2(V2l); 01497 ex V3m = lst_to_matrix2(V3l); 01498 01499 int k; 01500 for (unsigned int i=1; i < d; i++) 01501 { 01502 for (unsigned int j=1; j < d; j++) 01503 { 01504 if ( int(d) - int(i) - int(j) >= 1 ) 01505 { 01506 k = d - i - j; 01507 ex sum = (k*V1m + j*V2m + i*V3m)/d; 01508 ret.append(matrix_to_lst2(sum.evalm())); 01509 } 01510 } 01511 } 01512 // FIXME how should these be sorted ????? 01513 // ret = ret.sort(); 01514 return ret; 01515 } 01516 01517 lst bezier_ordinates(Line& line, unsigned int d) 01518 { 01519 01520 lst ret; 01521 ex V1 = line.vertex(0); 01522 ex V2 = line.vertex(1); 01523 01524 if (!GiNaC::is_a<lst>(V1)) 01525 { 01526 int k; 01527 for (unsigned int i=0; i <= d; i++) 01528 { 01529 k = d - i; 01530 ex sum = (k*V1 + i*V2)/d; 01531 ret.append(sum); 01532 } 01533 } 01534 else 01535 { 01536 01537 //FIXME: ugly conversion to matrix 01538 01539 lst V1l = ex_to<lst>(V1); 01540 lst V2l = ex_to<lst>(V2); 01541 01542 ex V1m = lst_to_matrix2(V1l); 01543 ex V2m = lst_to_matrix2(V2l); 01544 01545 int k; 01546 for (unsigned int i=0; i <= d; i++) 01547 { 01548 k = d - i; 01549 ex sum = (k*V1m + i*V2m)/d; 01550 ret.append(matrix_to_lst2(sum.evalm())); 01551 } 01552 // FIXME how should these be sorted ????? 01553 // ret = ret.sort(); 01554 } 01555 return ret; 01556 } 01557 01558 lst interior_coordinates(Line& line, unsigned int d) 01559 { 01560 01561 //FIXME: ugly conversion to matrix 01562 d = d+2; 01563 01564 lst ret; 01565 ex V1 = line.vertex(0); 01566 ex V2 = line.vertex(1); 01567 01568 lst V1l = ex_to<lst>(V1); 01569 lst V2l = ex_to<lst>(V2); 01570 01571 ex V1m = lst_to_matrix2(V1l); 01572 ex V2m = lst_to_matrix2(V2l); 01573 01574 int k; 01575 for (unsigned int i=1; i < d; i++) 01576 { 01577 k = d - i; 01578 ex sum = (k*V1m + i*V2m)/d; 01579 ret.append(matrix_to_lst2(sum.evalm())); 01580 } 01581 // FIXME how should these be sorted ????? 01582 // ret = ret.sort(); 01583 return ret; 01584 } 01585 01586 // FIXME barycenter_line, barycenter_triangle, barycenter_tetrahedron 01587 // all needs to determine which equations to use in a proper way 01588 01589 ex barycenter_line(ex p0, ex p1) 01590 { 01591 ex sol; 01592 01593 // 1D 01594 if (!GiNaC::is_a<lst>(p0)) 01595 { 01596 GiNaC::symbol b0("b0"), b1("b1"); 01597 ex eq1 = x == b0*p0 + b1*p1; 01598 ex eq2 = 1 == b0 + b1; 01599 sol = lsolve(lst(eq1, eq2), lst(b0, b1)); 01600 } 01601 else if (p0.nops() == 1 && p1.nops() == 1) 01602 { 01603 GiNaC::symbol b0("b0"), b1("b1"); 01604 ex eq1 = x == b0*p0.op(0) + b1*p1.op(0); 01605 ex eq2 = 1 == b0 + b1; 01606 sol = lsolve(lst(eq1, eq2), lst(b0, b1)); 01607 if ( sol == 0 ) 01608 { 01609 ex eq1 = y == b0*p0.op(1) + b1*p1.op(1); 01610 sol = lsolve(lst(eq1, eq2), lst(b0, b1)); 01611 } 01612 if ( sol == 0 ) 01613 { 01614 ex eq1 = z == b0*p0.op(2) + b1*p1.op(2); 01615 sol = lsolve(lst(eq1, eq2), lst(b0, b1)); 01616 } 01617 } 01618 //2D 01619 else if ( p0.nops() == 2 && p1.nops() == 2 ) 01620 { 01621 GiNaC::symbol b0("b0"), b1("b1"); 01622 ex eq1 = x == b0*p0.op(0) + b1*p1.op(0); 01623 ex eq3 = 1 == b0 + b1; 01624 sol = lsolve(lst(eq1, eq3), lst(b0, b1)); 01625 if (sol.nops() == 0) 01626 { 01627 ex eq2 = y == b0*p0.op(1) + b1*p1.op(1); 01628 sol = lsolve(lst(eq2, eq3), lst(b0, b1)); 01629 } 01630 } 01631 //3D 01632 else if ( p0.nops() == 3 && p1.nops() == 3 ) 01633 { 01634 GiNaC::symbol b0("b0"), b1("b1"); 01635 ex eq1 = x == b0*p0.op(0) + b1*p1.op(0); 01636 ex eq4 = 1 == b0 + b1; 01637 sol = lsolve(lst(eq1, eq4), lst(b0, b1)); 01638 if (sol.nops() == 0) 01639 { 01640 ex eq2 = y == b0*p0.op(1) + b1*p1.op(1); 01641 sol = lsolve(lst(eq2, eq4), lst(b0, b1)); 01642 } 01643 if (sol.nops() == 0) 01644 { 01645 ex eq3 = z == b0*p0.op(2) + b1*p1.op(2); 01646 sol = lsolve(lst(eq3, eq4), lst(b0, b1)); 01647 } 01648 } 01649 else 01650 { 01651 throw std::runtime_error("Could not compute the barycentric coordinates. Check the coordinates."); 01652 } 01653 01654 return sol; 01655 } 01656 01657 ex barycenter_triangle(ex p0, ex p1, ex p2) 01658 { 01659 ex sol; 01660 01661 // 2D 01662 if ( p0.nops() == 2 && p1.nops() == 2 && p2.nops() == 2) 01663 { 01664 GiNaC::symbol b0("b0"), b1("b1"), b2("b2"); 01665 ex eq1 = x == b0*p0.op(0) + b1*p1.op(0) + b2*p2.op(0); 01666 ex eq2 = y == b0*p0.op(1) + b1*p1.op(1) + b2*p2.op(1); 01667 ex eq3 = 1 == b0 + b1 + b2; 01668 01669 sol = lsolve(lst(eq1, eq2, eq3), lst(b0, b1, b2)); 01670 } 01671 // 3D 01672 else if ( p0.nops() == 3 && p1.nops() == 3 && p2.nops() == 3) 01673 { 01674 lst n1(p1.op(0) - p0.op(0), p1.op(1) - p0.op(1), p1.op(2) - p0.op(2)); 01675 lst n2 = lst(p2.op(0) - p0.op(0), p2.op(1) - p0.op(1), p2.op(2) - p0.op(2)); 01676 lst n = cross(n1,n2); 01677 01678 lst midpoint = lst((p0.op(0) + p1.op(0) + p2.op(0))/3, 01679 (p0.op(1) + p1.op(1) + p2.op(1))/3, 01680 (p0.op(2) + p1.op(2) + p2.op(2))/3); 01681 01682 ex p3 = lst(midpoint.op(0) + n.op(0), 01683 midpoint.op(1) + n.op(1), 01684 midpoint.op(2) + n.op(2)); 01685 01686 ex s = barycenter_tetrahedron(p0, p1, p2, p3); 01687 lst solution; 01688 for (unsigned int i=0; i<s.nops(); i++) 01689 { 01690 ex d = s.op(i).subs(x == p3.op(0)).subs(y == p3.op(1)).subs(z == p3.op(2)); 01691 d = d.rhs(); 01692 if ( GiNaC::is_a<GiNaC::numeric>(d)) 01693 { 01694 // FIXME: bad test, should use the toleranse variable set by CLN or something 01695 if ( GiNaC::abs(GiNaC::ex_to<GiNaC::numeric>(d)) < 10e-8) 01696 { 01697 solution.append(s.op(i)); 01698 } 01699 } 01700 else 01701 { 01702 if ( d.is_zero() ) 01703 { 01704 solution.append(s.op(i)); 01705 } 01706 } 01707 } 01708 sol = solution; 01709 } 01710 else 01711 { 01712 throw std::runtime_error("Could not compute the barycentric coordinates. Check the coordinates."); 01713 } 01714 01715 return sol; 01716 } 01717 01718 ex barycenter_tetrahedron(ex p0, ex p1, ex p2, ex p3) 01719 { 01720 GiNaC::symbol b0("b0"), b1("b1"), b2("b2"), b3("b3"); 01721 01722 // 3D 01723 ex eq1 = x == b0*p0.op(0) + b1*p1.op(0) + b2*p2.op(0) + b3*p3.op(0); 01724 ex eq2 = y == b0*p0.op(1) + b1*p1.op(1) + b2*p2.op(1) + b3*p3.op(1); 01725 ex eq3 = z == b0*p0.op(2) + b1*p1.op(2) + b2*p2.op(2) + b3*p3.op(2); 01726 ex eq4 = 1 == b0 + b1 + b2 +b3; 01727 01728 ex sol = lsolve(lst(eq1, eq2, eq3, eq4), lst(b0, b1, b2, b3)); 01729 01730 return sol; 01731 01732 } 01733 01734 ex barycenter(Simplex& simplex) 01735 { 01736 if (nsd != simplex.no_vertices()-1) 01737 { 01738 throw std::runtime_error("Could not compute the barycentric coordinates. Not implemented yet for simplices with no_vertices != nsd +1."); 01739 } 01740 01741 // put symbols in lst 01742 ex b = get_symbolic_vector(simplex.no_vertices(), "b"); 01743 lst symbols; 01744 for (unsigned int i=0; i<b.nops(); i++) 01745 { 01746 symbols.append(b.op(i)); 01747 } 01748 01749 // put equations in lst 01750 lst eqs; 01751 for (unsigned int i=0; i<nsd; i++) 01752 { 01753 ex sum = 0; 01754 for (unsigned int k=0; k< simplex.no_vertices(); k++) 01755 { 01756 sum += b.op(k)*simplex.vertex(k).op(i); 01757 } 01758 ex eqi = p[i] == sum; 01759 eqs.append(eqi); 01760 } 01761 01762 // last eq, sum = 1 01763 ex sum = 0; 01764 for (unsigned int i=0; i<symbols.nops(); i++) 01765 { 01766 sum += symbols.op(i); 01767 } 01768 ex last_eq = 1 == sum; 01769 eqs.append(last_eq); 01770 01771 // solve equations 01772 ex sol = lsolve(eqs, symbols); 01773 return sol; 01774 } 01775 01776 ex bernstein(unsigned int order, Polygon& p, const string & a) 01777 { 01778 01779 if ( order < 0 ) 01780 { 01781 throw(std::logic_error("Can not create polynomials of order less than 0!")); 01782 } 01783 01784 ex ret; // ex to return 01785 int dof; // degrees of freedom 01786 ex A; // ex holding the coefficients a_0 .. a_dof 01787 lst basis; 01788 01789 if ( p.str().find("Line") != string::npos ) 01790 { 01791 ex bary = barycenter_line(p.vertex(0), p.vertex(1)); 01792 ex b0= bary.op(0).rhs(); 01793 ex b1= bary.op(1).rhs(); 01794 dof = order+1; 01795 A = get_symbolic_matrix(1,dof, a); 01796 int o=0; 01797 for (GiNaC::const_iterator i = A.begin(); i != A.end(); ++i) 01798 { 01799 ex scale = GiNaC::binomial(order,o); 01800 ret += (*i)*scale*pow(b0,o)*pow(b1,order-o); 01801 basis.append(scale*pow(b0,o)*pow(b1,order-o)); 01802 o++; 01803 } 01804 } 01805 else if ( p.str().find("Triangle") != string::npos ) 01806 { 01807 01808 dof = (order+1)*(order+2)/2; 01809 A = get_symbolic_matrix(1, dof , a); 01810 01811 ex bary = barycenter_triangle(p.vertex(0), p.vertex(1), p.vertex(2)); 01812 ex b0= bary.op(0).rhs(); 01813 ex b1= bary.op(1).rhs(); 01814 ex b2= bary.op(2).rhs(); 01815 01816 size_t i=0; 01817 for (unsigned int o1 = 0; o1 <= order; o1++) 01818 { 01819 for (unsigned int o2 = 0; o2 <= order; o2++) 01820 { 01821 for (unsigned int o3 = 0; o3 <= order; o3++) 01822 { 01823 if ( o1 + o2 + o3 == order ) 01824 { 01825 ex scale = (GiNaC::factorial(order)/(GiNaC::factorial(o1)*GiNaC::factorial(o2)*GiNaC::factorial(o3))); 01826 ret += A.op(i)*scale*pow(b0,o1)*pow(b1,o2)*pow(b2,o3); 01827 01828 basis.append(scale*pow(b0,o1)*pow(b1,o2)*pow(b2,o3)); 01829 i++; 01830 } 01831 } 01832 } 01833 } 01834 } 01835 01836 else if ( p.str().find("Tetrahedron") != string::npos ) 01837 { 01838 01839 dof = 0; 01840 for (unsigned int j=0; j<= order; j++) 01841 { 01842 dof += (j+1)*(j+2)/2; 01843 } 01844 A = get_symbolic_matrix(1, dof , a); 01845 01846 ex bary = barycenter_tetrahedron(p.vertex(0), p.vertex(1), p.vertex(2), p.vertex(3)); 01847 ex b0= bary.op(0).rhs(); 01848 ex b1= bary.op(1).rhs(); 01849 ex b2= bary.op(2).rhs(); 01850 ex b3= bary.op(3).rhs(); 01851 01852 size_t i=0; 01853 for (unsigned int o1 = 0; o1 <= order; o1++) 01854 { 01855 for (unsigned int o2 = 0; o2 <= order; o2++) 01856 { 01857 for (unsigned int o3 = 0; o3 <= order; o3++) 01858 { 01859 for (unsigned int o4 = 0; o4 <= order; o4++) 01860 { 01861 if ( o1 + o2 + o3 + o4 == order ) 01862 { 01863 ex scale = (GiNaC::factorial(order)/(GiNaC::factorial(o1)*GiNaC::factorial(o2)*GiNaC::factorial(o3)*GiNaC::factorial(o4))); 01864 ret += A.op(i)*scale*pow(b0,o1)*pow(b1,o2)*pow(b2,o3)*pow(b3,o4); 01865 basis.append(scale*pow(b0,o1)*pow(b1,o2)*pow(b2,o3)*pow(b3,o4)); 01866 i++; 01867 } 01868 } 01869 } 01870 } 01871 } 01872 } 01873 01874 else if (p.str() == "Simplex" || p.str() == "ReferenceSimplex") 01875 { 01876 01877 throw std::runtime_error("Not implemented yet."); 01878 // ex bary = barycenter(p); 01879 } 01880 return lst(ret,matrix_to_lst2(A),basis); 01881 } 01882 01883 lst bernsteinv(unsigned int no_fields, unsigned int order, Polygon& p, const string & a) 01884 { 01885 01886 if ( order < 0 ) 01887 { 01888 throw(std::logic_error("Can not create polynomials of order less than 0!")); 01889 } 01890 01891 lst ret1; // contains the polynom 01892 lst ret2; // contains the coefficients 01893 lst ret3; // constains the basis functions 01894 lst basis_tmp; 01895 for (unsigned int i=0; i< no_fields; i++) 01896 { 01897 lst basis; 01898 std::ostringstream s; 01899 s <<a<<""<<i<<"_"; 01900 ex pol = bernstein(order, p, s.str()); 01901 ret1.append(pol.op(0)); 01902 ret2.append(pol.op(1)); 01903 basis_tmp = ex_to<lst>(pol.op(2)); 01904 for (lst::const_iterator basis_iterator = basis_tmp.begin(); 01905 basis_iterator != basis_tmp.end(); ++basis_iterator) 01906 { 01907 lst tmp_lst; 01908 for (unsigned int d=1; d<=no_fields; d++) tmp_lst.append(0); 01909 tmp_lst.let_op(i) = (*basis_iterator); 01910 ret3.append(tmp_lst); 01911 } 01912 } 01913 return lst(ret1,ret2,ret3); 01914 01915 } 01916 01917 lst normal(Tetrahedron& tetrahedron, unsigned int i) 01918 { 01919 // Normal as defined by Maries note 01920 Triangle triangle = tetrahedron.triangle(i); 01921 lst vertex_i = ex_to<lst>(tetrahedron.vertex(i)); 01922 lst vertex_0 = ex_to<lst>(triangle.vertex(0)); 01923 lst vertex_1 = ex_to<lst>(triangle.vertex(1)); 01924 lst vertex_2 = ex_to<lst>(triangle.vertex(2)); 01925 01926 lst n1(vertex_1.op(0) - vertex_0.op(0), 01927 vertex_1.op(1) - vertex_0.op(1), 01928 vertex_1.op(2) - vertex_0.op(2)); 01929 01930 lst n2(vertex_2.op(0) - vertex_0.op(0), 01931 vertex_2.op(1) - vertex_0.op(1), 01932 vertex_2.op(2) - vertex_0.op(2)); 01933 01934 /* 01935 lst n3(vertex_0.op(0) - vertex_i.op(0), 01936 vertex_0.op(1) - vertex_i.op(1), 01937 vertex_0.op(2) - vertex_i.op(2)); 01938 */ 01939 01940 lst n4 = cross(n1,n2); 01941 /* 01942 ex nn = inner(n3, n4); 01943 int sign = 1; 01944 if ( is_a<numeric>(nn)) { 01945 if ( nn > 0 ) { 01946 sign = 1; 01947 } else if ( nn < 0) { 01948 sign = -1; 01949 } else { 01950 sign = 0; 01951 } 01952 } 01953 */ 01954 01955 ex norm = sqrt(pow(n4.op(0),2) 01956 + pow(n4.op(1),2) 01957 + pow(n4.op(2),2)); 01958 01959 n4.let_op(0) = n4.op(0)/norm; 01960 n4.let_op(1) = n4.op(1)/norm; 01961 n4.let_op(2) = n4.op(2)/norm; 01962 01963 return n4; 01964 01965 } 01966 01967 lst normal(Triangle& triangle, unsigned int i) 01968 { 01969 Line line = triangle.line(i); 01970 lst vertex_i = ex_to<lst>(triangle.vertex(i)); 01971 lst vertex_0 = ex_to<lst>(line.vertex(0)); 01972 lst vertex_1 = ex_to<lst>(line.vertex(1)); 01973 01974 /* 01975 lst n1 = lst (- (vertex_1.op(1) - vertex_0.op(1)), vertex_1.op(0) - vertex_0.op(0) ); 01976 lst n2 = lst (vertex_0.op(0) - vertex_i.op(0), vertex_0.op(1) - vertex_i.op(1)); 01977 01978 ex nn = inner(n1, n2); 01979 int sign = 1; 01980 / * 01981 if ( is_a<numeric>(nn)) { 01982 if ( nn > 0 ) { 01983 sign = 1; 01984 } else if ( nn < 0) { 01985 sign = -1; 01986 } else { 01987 sign = 0; 01988 } 01989 } 01990 01991 ex norm = sqrt(pow(n1.op(0),2) + pow(n1.op(1),2)); 01992 n1.let_op(0) = sign*n1.op(0)/norm; 01993 n1.let_op(1) = sign*n1.op(1)/norm; 01994 */ 01995 01996 // normal vector as Marie has defined them 01997 lst n1 = lst ( (vertex_1.op(1) - vertex_0.op(1)), 01998 -(vertex_1.op(0) - vertex_0.op(0)) ); 01999 02000 ex norm = sqrt(pow(n1.op(0),2) + pow(n1.op(1),2)); 02001 n1.let_op(0) = n1.op(0)/norm; 02002 n1.let_op(1) = n1.op(1)/norm; 02003 02004 return n1; 02005 } 02006 02007 lst tangent(Triangle& triangle, unsigned int i) 02008 { 02009 /* 02010 Line line = triangle.line(i); 02011 //FIXME: 5 lines to compute the tangent vector, these should 02012 // be put somewhere else. 02013 GiNaC::symbol t("t"); 02014 ex line_repr = line.repr(t); 02015 ex t1 = line_repr.op(0).rhs().coeff(t,1); 02016 ex t2 = line_repr.op(1).rhs().coeff(t,1); 02017 ex norm = sqrt(pow(t1,2) + pow(t2,2)); 02018 lst tangent = lst(t1/norm,t2/norm); 02019 return tangent; 02020 */ 02021 /* 02022 ex t1, t2; 02023 if ( i == 0 ) { 02024 t1 = triangle.vertex(2).op(0) - triangle.vertex(1).op(0); 02025 t2 = triangle.vertex(2).op(1) - triangle.vertex(1).op(1); 02026 } else if ( i == 1 ) { 02027 t1 = triangle.vertex(0).op(0) - triangle.vertex(2).op(0); 02028 t2 = triangle.vertex(0).op(1) - triangle.vertex(2).op(1); 02029 } else if ( i == 2 ) { 02030 t1 = triangle.vertex(1).op(0) - triangle.vertex(0).op(0); 02031 t2 = triangle.vertex(1).op(1) - triangle.vertex(0).op(1); 02032 } else { 02033 throw(std::out_of_range("The side index is out of range!")); 02034 } 02035 */ 02036 Line line = triangle.line(i); 02037 ex t1 = line.vertex(1).op(0) - line.vertex(0).op(0); 02038 ex t2 = line.vertex(1).op(1) - line.vertex(0).op(1); 02039 02040 ex norm = sqrt(pow(t1,2) + pow(t2,2)); 02041 lst tangent = lst(t1/norm,t2/norm); 02042 return tangent; 02043 02044 } 02045 02046 } //namespace SyFi