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