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