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