SyFi 0.3
Lagrange.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 "Lagrange.h"
00007 #include "ElementComputations.h"
00008 #include "tools.h"
00009 
00010 using std::cout;
00011 using std::endl;
00012 using std::string;
00013 
00014 namespace SyFi
00015 {
00016 
00017         Lagrange::Lagrange() : StandardFE()
00018         {
00019                 description = "Lagrange";
00020         }
00021 
00022         Lagrange::Lagrange(Polygon& p, unsigned int order) : StandardFE(p, order)
00023         {
00024                 compute_basis_functions();
00025         }
00026 
00027         void Lagrange:: compute_basis_functions()
00028         {
00029 
00030                 // NOTE: in the below code dof(i) is not used to
00031                 // determine the basis functions
00032 
00033                 // remove previously computed basis functions and dofs
00034                 Ns.clear();
00035                 dofs.clear();
00036 
00037                 if ( order < 1 )
00038                 {
00039                         throw(std::logic_error("Lagrangian elements must be of order 1 or higher."));
00040                 }
00041 
00042                 if ( p == NULL )
00043                 {
00044                         throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00045                 }
00046 
00047                 GiNaC::lst equations;
00048                 GiNaC::lst variables;
00049                 GiNaC::ex polynom;
00050 
00051                 if (p->str().find("ReferenceLine") != string::npos)
00052                 {
00053 
00054                         description = istr("Lagrange_", order) + "_1D";
00055 
00056                         // Look at the case with the Triangle for a documented code
00057 
00058                         //    polynom = pol(order, 1, "a");
00059                         //    variables = coeffs(polynom);
00060                         GiNaC::ex polynom_space = bernstein(order, *p, "a");
00061                         //    GiNaC::ex polynom_space = pol(order, 1, "a");
00062                         polynom = polynom_space.op(0);
00063                         variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00064 
00065                         GiNaC::ex increment = GiNaC::numeric(1,order);
00066 
00067                         GiNaC::ex Nj;
00068                         for (GiNaC::ex p=0; p<= 1 ; p += increment )
00069                         {
00070                                 GiNaC::ex eq = polynom == GiNaC::numeric(0);
00071                                 equations.append(eq.subs(x == p));
00072                                 dofs.insert(dofs.end(), GiNaC::lst(p));
00073                         }
00074                 }
00075                 else if (p->str().find("Line") != string::npos  )
00076                 {
00077                         description = istr("Lagrange_", order) + "_1D";
00078 
00079                         // Look at the case with the Triangle for a documented code
00080 
00081                         //    polynom = pol(order, 1, "a");
00082                         //    variables = coeffs(polynom);
00083                         GiNaC::ex polynom_space = bernstein(order, *p, "a");
00084                         //    GiNaC::ex polynom_space = pol(order, 1, "a");
00085                         polynom = polynom_space.op(0);
00086                         variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00087 
00088                         Polygon& pp = *p;
00089                         Line& l = (Line&) pp;
00090                         GiNaC::lst points = bezier_ordinates(l,order);
00091 
00092                         GiNaC::ex Nj;
00093                         for (unsigned int i=1; i<= points.nops() ; i++ )
00094                         {
00095                                 GiNaC::ex point = points.op(i-1);
00096                                 GiNaC::ex eq = polynom == GiNaC::numeric(0);
00097                                 if (point.nops() == 0) eq = eq.subs(x == point);
00098                                 if (point.nops() > 0) eq = eq.subs(x == point.op(0));
00099                                 if (point.nops() > 1) eq = eq.subs(y == point.op(1));
00100                                 if (point.nops() > 2) eq = eq.subs(z == point.op(2));
00101                                 equations.append(eq);
00102                                 dofs.insert(dofs.end(), GiNaC::lst(point));
00103                         }
00104 
00105                 }
00106                 else if (p->str().find("ReferenceTriangle") != string::npos  )
00107                 {
00108 
00109                         description = istr("Lagrange_", order) + "_2D";
00110 
00111                         // Look at the case with the Triangle for a documented code
00112 
00113                         //    polynom = pol(order, 2, "b");
00114                         //    variables = coeffs(polynom);
00115                         GiNaC::ex polynom_space = bernstein(order, *p, "b");
00116                         //    GiNaC::ex polynom_space = pol(order, 2, "a");
00117                         polynom = polynom_space.op(0);
00118                         variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00119 
00120                         GiNaC::ex increment = GiNaC::numeric(1,order);
00121 
00122                         GiNaC::ex Nj;
00123                         GiNaC::numeric one = 1;
00124                         for (GiNaC::ex q=0; q<= one  ; q += increment )
00125                         {
00126                                 for (GiNaC::ex p=0; p<= one-q ; p += increment )
00127                                 {
00128                                         GiNaC::ex eq = polynom == GiNaC::numeric(0);
00129                                         equations.append(eq.subs(GiNaC::lst(x == p, y == q)));
00130                                         dofs.insert(dofs.end(), GiNaC::lst(p,q));
00131                                 }
00132                         }
00133                 }
00134                 else if ( p->str().find("Triangle") != string::npos)
00135                 {
00136 
00137                         description = istr("Lagrange_", order) + "_2D";
00138 
00139                         // Look HERE for the documented code
00140 
00141                         //        GiNaC::ex polynom_space = pol(order, 2, "a");
00142                         GiNaC::ex polynom_space = bernstein(order, *p, "b");
00143                         // the polynomial spaces on the form:
00144                         // first item:     a0 + a1*x + a2*y + a3*x^2 + a4*x*y ...     the polynom
00145                         // second item:    a0, a1, a2, ...                            the coefficents
00146                         // third  item     1, x, y, x^2, ..                           the basis
00147                         // Could also do:
00148                         // GiNaC::ex polynom_space = bernstein(order, t, "a");
00149 
00150                         polynom = polynom_space.op(0);
00151                         variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00152 
00153                         GiNaC::ex Nj;
00154                         Polygon& pp = *p;
00155                         Triangle& t = (Triangle&) pp;
00156                         // The bezier ordinates (in which the basis function should be either 0 or 1)
00157                         GiNaC::lst points = bezier_ordinates(t,order);
00158 
00159                         for (unsigned int i=1; i<= points.nops() ; i++ )
00160                         {
00161                                 GiNaC::ex point = points.op(i-1);
00162                                 GiNaC::ex eq = polynom == GiNaC::numeric(0);
00163                                 equations.append(eq.subs(GiNaC::lst(x == point.op(0) , y == point.op(1))));
00164                                 dofs.insert(dofs.end(), GiNaC::lst(point.op(0),point.op(1)));
00165                         }
00166                 }
00167 
00168                 else if ( p->str().find("ReferenceRectangle") != string::npos)
00169                 {
00170 
00171                         description = istr("Lagrange_", order) + "_2D";
00172 
00173                         // create 1D element, then create tensor product
00174                         ReferenceLine line;
00175                         Lagrange fe(line, order);
00176 
00177                         for (unsigned int i=0; i< fe.nbf(); i++)
00178                         {
00179                                 for (unsigned int j=0; j< fe.nbf(); j++)
00180                                 {
00181                                         Ns.insert(Ns.end(), fe.N(i)*fe.N(j).subs(x==y));
00182                                         dofs.insert(dofs.end(), GiNaC::lst(fe.dof(i).op(0), fe.dof(j).op(0)));
00183                                 }
00184                         }
00185                         return;
00186 
00187                         /* OLD CODE
00188 
00189                            GiNaC::ex polynom_space = legendre(order, 2, "b");
00190 
00191                            polynom = polynom_space.op(0);
00192                            variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00193 
00194                         GiNaC::ex Nj;
00195                         Polygon& pp = *p;
00196                         ReferenceRectangle& b = (ReferenceRectangle&) pp;
00197 
00198                         int no_points = (order+1)*(order+1);
00199                         GiNaC::ex increment = GiNaC::numeric(1,order);
00200 
00201                         GiNaC::numeric one=1.0;
00202 
00203                         for (GiNaC::ex q=0; q <= one  ; q += increment )
00204                         {
00205                         for (GiNaC::ex p=0; p <= one ; p += increment )
00206                         {
00207                         GiNaC::ex eq = polynom == GiNaC::numeric(0);
00208                         equations.append(eq.subs(GiNaC::lst(x == p, y == q)));
00209                         dofs.push_back(GiNaC::lst(p,q));
00210                         }
00211                         }
00212                         */
00213                 }
00214                 else if ( p->str().find("ReferenceTetrahedron") != string::npos)
00215                 {
00216 
00217                         description = istr("Lagrange_", order) + "_3D";
00218 
00219                         // Look at the case with the Triangle for a documented code
00220 
00221                         //    polynom = pol(order, 3, "b");
00222                         //    GiNaC::ex polynom_space = pol(order, 3, "a");
00223                         GiNaC::ex polynom_space = bernstein(order, *p, "b");
00224                         polynom = polynom_space.op(0);
00225                         variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00226 
00227                         int nno =0;
00228                         for (unsigned int j=0; j<= order; j++)
00229                         {
00230                                 nno += (j+1)*(j+2)/2;
00231                         }
00232 
00233                         GiNaC::ex increment = GiNaC::numeric(1,order);
00234 
00235                         GiNaC::ex Nj;
00236                         for (GiNaC::ex r=0; r<= 1 ; r += increment )
00237                         {
00238                                 for (GiNaC::ex q=0; q<= 1-r ; q += increment )
00239                                 {
00240                                         for (GiNaC::ex p=0; p<= 1-r-q ; p += increment )
00241                                         {
00242                                                 GiNaC::ex eq = polynom == GiNaC::numeric(0);
00243                                                 equations.append(eq.subs(GiNaC::lst(x == p, y == q, z == r )));
00244                                                 dofs.push_back(GiNaC::lst(p,q,r));
00245                                         }
00246                                 }
00247                         }
00248                 }
00249                 else if ( p->str().find("Tetrahedron") != string::npos)
00250                 {
00251 
00252                         description = istr("Lagrange_", order) + "_3D";
00253 
00254                         // Look at the case with the Triangle for a documented code
00255 
00256                         GiNaC::ex polynom_space = pol(order, 3, "a");
00257                         //    GiNaC::ex polynom_space = bernstein(order, *p, "b");
00258                         polynom = polynom_space.op(0);
00259                         variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00260 
00261                         GiNaC::ex increment = GiNaC::numeric(1,order);
00262 
00263                         GiNaC::ex Nj;
00264                         Polygon& pp = *p;
00265                         Tetrahedron& t = (Tetrahedron&) pp;
00266                         GiNaC::lst points = bezier_ordinates(t,order);
00267                         for (unsigned int i=1; i<= points.nops() ; i++ )
00268                         {
00269                                 GiNaC::ex point = points.op(i-1);
00270                                 GiNaC::ex eq = polynom == GiNaC::numeric(0);
00271                                 equations.append(eq.subs(GiNaC::lst(x == point.op(0) , y == point.op(1), z == point.op(2))));
00272                                 dofs.push_back(GiNaC::lst(point.op(0),point.op(1),point.op(2)));
00273                         }
00274                 }
00275                 else if ( p->str().find("ReferenceBox") != string::npos)
00276                 {
00277 
00278                         description = istr("Lagrange_", order) + "_3D";
00279 
00280                         ReferenceLine line;
00281                         Lagrange fe(line, order);
00282 
00283                         for (unsigned int i=0; i< fe.nbf(); i++)
00284                         {
00285                                 for (unsigned int j=0; j< fe.nbf(); j++)
00286                                 {
00287                                         for (unsigned int k=0; k< fe.nbf(); k++)
00288                                         {
00289                                                 Ns.insert(Ns.end(), fe.N(i)*fe.N(j).subs(x==y)*fe.N(k).subs(x==z));
00290                                                 dofs.insert(dofs.end(), GiNaC::lst(fe.dof(i).op(0), fe.dof(j).op(0), fe.dof(k).op(0)));
00291                                         }
00292                                 }
00293                         }
00294                         return;
00295 
00296                         /* OLD CODE
00297 
00298                            GiNaC::ex polynom_space = legendre(order, 3, "b");
00299 
00300                            polynom = polynom_space.op(0);
00301                            variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00302 
00303                         GiNaC::ex Nj;
00304                         Polygon& pp = *p;
00305                         ReferenceRectangle& b = (ReferenceRectangle&) pp;
00306 
00307                         int no_points = (order+1)*(order+1)*(order+1);
00308                         GiNaC::ex increment = GiNaC::numeric(1,order);
00309 
00310                         GiNaC::numeric one = 1;
00311                         for (GiNaC::ex p=0; p <= one ; p += increment) {
00312                         for (GiNaC::ex q=0; q <= one  ; q += increment) {
00313                         for (GiNaC::ex r=0; r <= one ; r += increment) {
00314                         GiNaC::ex eq = polynom == GiNaC::numeric(0);
00315                         equations.append(eq.subs(GiNaC::lst(x == p, y == q, z==r)));
00316                         dofs.push_back(GiNaC::lst(p,q,r));
00317                         }
00318                         }
00319                         }
00320 
00321                         / *
00322                                 GiNaC::ex subs = lsolve(equations, variables);
00323                                 Nj = polynom.subs(subs);
00324                                 Ns.push_back(Nj);
00325                                 */
00326                 }
00327 
00328                 // invert the matrix:
00329                 // GiNaC has a bit strange way to invert a matrix.
00330                 // It solves the system AA^{-1} = Id.
00331                 // It seems that this way is the only way to do
00332                 // properly with the solve_algo::gauss flag.
00333                 //
00334 
00335                 //    std::cout <<"no variables "<<variables.nops()<<std::endl;
00336                 //    std::cout <<"no equations "<<equations.nops()<<std::endl;
00337                 //    print(equations);
00338                 //    print(variables);
00339                 GiNaC::matrix b; GiNaC::matrix A;
00340                 matrix_from_equations(equations, variables, A, b);
00341 
00342                 unsigned int ncols = A.cols();
00343                 GiNaC::matrix vars_sq(ncols, ncols);
00344 
00345                 // matrix of symbols
00346                 for (unsigned r=0; r<ncols; ++r)
00347                         for (unsigned c=0; c<ncols; ++c)
00348                                 vars_sq(r, c) = GiNaC::symbol();
00349 
00350                 GiNaC::matrix id(ncols, ncols);
00351 
00352                 // identity
00353                 const GiNaC::ex _ex1(1);
00354                 for (unsigned i=0; i<ncols; ++i)
00355                         id(i, i) = _ex1;
00356 
00357                 // invert the matrix
00358                 GiNaC::matrix m_inv(ncols, ncols);
00359                 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss);
00360 
00361                 for (unsigned int i=0; i<dofs.size(); i++)
00362                 {
00363                         b.let_op(i) = GiNaC::numeric(1);
00364                         GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b));
00365 
00366                         GiNaC::lst subs;
00367                         for (unsigned int ii=0; ii<xx.nops(); ii++)
00368                         {
00369                                 subs.append(variables.op(ii) == xx.op(ii));
00370                         }
00371                         GiNaC::ex Nj = polynom.subs(subs);
00372                         Ns.insert(Ns.end(), Nj);
00373                         b.let_op(i) = GiNaC::numeric(0);
00374                 }
00375         }
00376 
00377         // ------------VectorLagrange ---
00378 
00379         VectorLagrange::VectorLagrange() : StandardFE()
00380         {
00381                 description = "VectorLagrange";
00382         }
00383 
00384         VectorLagrange::VectorLagrange(Polygon& p, unsigned int order, unsigned int size_) : StandardFE(p, order)
00385         {
00386                 size = size_ < 0 ? nsd: size_;
00387                 compute_basis_functions();
00388         }
00389 
00390         void VectorLagrange:: compute_basis_functions()
00391         {
00392 
00393                 // remove previously computed basis functions and dofs
00394                 Ns.clear();
00395                 dofs.clear();
00396 
00397                 if ( order < 1 )
00398                 {
00399                         throw(std::logic_error("Lagrangian elements must be of order 1 or higher."));
00400                 }
00401 
00402                 if ( p == NULL )
00403                 {
00404                         throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00405                 }
00406 
00407                 if ( size == 0)
00408                 {
00409                         throw(std::logic_error("You need to set the size of the vector before the basisfunctions can be computed"));
00410                 }
00411 
00412                 Lagrange fe;
00413                 fe.set_order(order);
00414                 fe.set_polygon(*p);
00415                 fe.compute_basis_functions();
00416                 GiNaC::lst zero_list;
00417                 for (unsigned int s=1; s<= size ; s++)
00418                 {
00419                         zero_list.append(0);
00420                 }
00421 
00422                 for (unsigned int s=0; s< size ; s++)
00423                 {
00424                         for (unsigned int i=0; i< fe.nbf() ; i++)
00425                         {
00426                                 GiNaC::lst Nis = zero_list;
00427                                 Nis.let_op(s) = fe.N(i);
00428                                 GiNaC::ex Nmat = GiNaC::matrix(size,1,Nis);
00429                                 Ns.insert(Ns.end(), Nmat);
00430 
00431                                 GiNaC::lst dof = GiNaC::lst(fe.dof(i), s) ;
00432                                 dofs.insert(dofs.end(), dof);
00433                         }
00434                 }
00435 
00436                 description = "Vector" + fe.str();
00437         }
00438 
00439         void VectorLagrange:: set_size(unsigned int size_)
00440         {
00441                 size = size_;
00442         }
00443 
00444         // ------------TensorLagrange ---
00445 
00446         TensorLagrange::TensorLagrange() : StandardFE()
00447         {
00448                 description = "TensorLagrange";
00449         }
00450 
00451         TensorLagrange::TensorLagrange(Polygon& p, unsigned int order, unsigned int size_) : StandardFE(p, order)
00452         {
00453                 size = size_ < 0 ? nsd: size_;
00454                 compute_basis_functions();
00455         }
00456 
00457         void TensorLagrange:: compute_basis_functions()
00458         {
00459 
00460                 // remove previously computed basis functions and dofs
00461                 Ns.clear();
00462                 dofs.clear();
00463 
00464                 if ( order < 1 )
00465                 {
00466                         throw(std::logic_error("Lagrangian elements must be of order 1 or higher."));
00467                 }
00468 
00469                 if ( p == NULL )
00470                 {
00471                         throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00472                 }
00473 
00474                 if ( size == 0)
00475                 {
00476                         throw(std::logic_error("You need to set the size of the vector before the basisfunctions can be computed"));
00477                 }
00478 
00479                 Lagrange fe;
00480                 fe.set_order(order);
00481                 fe.set_polygon(*p);
00482                 fe.compute_basis_functions();
00483                 GiNaC::lst zero_list;
00484                 for (unsigned int s=1; s<= size*size ; s++)
00485                 {
00486                         zero_list.append(0);
00487                 }
00488 
00489                 for (unsigned int r=0; r< size ; r++)
00490                 {
00491                         for (unsigned int s=0; s< size ; s++)
00492                         {
00493                                 for (unsigned int i=0; i< fe.nbf() ; i++)
00494                                 {
00495                                         GiNaC::lst Nis = zero_list;
00496                                         Nis.let_op((size)*r + s) = fe.N(i);
00497                                         GiNaC::ex Nmat = GiNaC::matrix(size,size,Nis);
00498                                         Ns.insert(Ns.end(), Nmat);
00499 
00500                                         GiNaC::lst dof = GiNaC::lst(fe.dof(i), r, s) ;
00501                                         dofs.insert(dofs.end(), dof);
00502                                 }
00503                         }
00504                 }
00505 
00506                 description = "Tensor" + fe.str();
00507         }
00508 
00509         void TensorLagrange:: set_size(unsigned int size_)
00510         {
00511                 size = size_;
00512         }
00513 
00514         GiNaC::ex lagrange(unsigned int order, Polygon& p, const std::string & a)
00515         {
00516                 if ( order < 1 )
00517                 {
00518                         throw(std::logic_error("Can not create polynomials of order less than 1!"));
00519                 }
00520 
00521                 GiNaC::ex A;
00522                 GiNaC::ex ret;
00523                 GiNaC::lst basis;
00524 
00525                 Lagrange fe(p,order);
00526                 A = get_symbolic_matrix(1, fe.nbf(), a);
00527 
00528                 for (unsigned int i=0; i<fe.nbf(); i++)
00529                 {
00530                         ret += A.op(i)*fe.N(i);
00531                         basis.append(fe.N(i));
00532                 }
00533                 return GiNaC::lst(ret,matrix_to_lst2(A),basis);
00534         }
00535 
00536         GiNaC::lst lagrangev(unsigned int no_fields, unsigned int order, Polygon& p, const std::string & a)
00537         {
00538                 if ( order < 1 )
00539                 {
00540                         throw(std::logic_error("Can not create polynomials of order less than 1!"));
00541                 }
00542 
00543                 GiNaC::ex A;
00544                 GiNaC::ex ret;
00545                 GiNaC::lst basis;
00546 
00547                 VectorLagrange fe(p,order);
00548                 A = get_symbolic_matrix(1, fe.nbf(), a);
00549 
00550                 for (unsigned int i=0; i<fe.nbf(); i++)
00551                 {
00552                         ret += A.op(i)*fe.N(i);
00553                         basis.append(fe.N(i));
00554                 }
00555                 return GiNaC::lst(ret,matrix_to_lst2(A),basis);
00556         }
00557 
00558 }                                                                // namespace SyFi
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines