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