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