SyFi 0.3
Robust.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 "Robust.h"
00005 #include <fstream>
00006 #include "tools.h"
00007 
00008 using std::cout;
00009 using std::endl;
00010 using std::string;
00011 using GiNaC::exmap;
00012 
00013 namespace SyFi
00014 {
00015 
00016         Robust:: Robust() : StandardFE()
00017         {
00018                 description = "Robust";
00019         }
00020 
00021         Robust:: Robust(Polygon& p, int unsigned order, bool pointwise_) : StandardFE(p, order)
00022         {
00023                 pointwise = pointwise_;
00024                 compute_basis_functions();
00025         }
00026 
00027         void Robust:: compute_basis_functions()
00028         {
00029 
00030                 // remove previously computed basis functions and dofs
00031                 Ns.clear();
00032                 dofs.clear();
00033 
00034                 GiNaC::ex nsymb = GiNaC::symbol("n");
00035                 GiNaC::ex tsymb = GiNaC::symbol("t");
00036 
00037                 if ( p == NULL )
00038                 {
00039                         throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00040                 }
00041 
00042                 if ( p->str().find("Line") != string::npos )
00043                 {
00044 
00045                         cout <<"Can not define the Robust element on a line"<<endl;
00046 
00047                 }
00048                 else if ( p->str().find("Triangle") != string::npos )
00049                 {
00050 
00051                         if (pointwise)
00052                         {
00053 
00054                                 description = "Robust_2D";
00055 
00056                                 Triangle& triangle = (Triangle&)(*p);
00057                                 GiNaC::lst equations;
00058                                 GiNaC::lst variables;
00059                                 GiNaC::ex V_space = bernsteinv(2, order+3, triangle, "a");
00060                                 GiNaC::ex V = V_space.op(0);
00061                                 GiNaC::ex V_vars = V_space.op(1);
00062 
00063                                 GiNaC::ex divV = div(V);
00064                                 exmap basis2coeff = pol2basisandcoeff(divV);
00065                                 exmap::iterator iter;
00066 
00067                                 // div constraints:
00068                                 for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
00069                                 {
00070                                         GiNaC::ex basis = (*iter).first;
00071                                         GiNaC::ex coeff= (*iter).second;
00072                                         if ( coeff != 0 )
00073                                         {
00074                                         }
00075                                         if ( coeff != 0 && ( basis.degree(x) + basis.degree(y)  > int(order) )  )
00076                                         {
00077                                                 equations.append( coeff == 0 );
00078                                         }
00079                                 }
00080 
00081                                 // normal constraints on edges:
00082                                 for (int i=0; i< 3; i++)
00083                                 {
00084                                         Line line = triangle.line(i);
00085                                         GiNaC::symbol s("s");
00086                                         GiNaC::lst normal_vec = normal(triangle, i);
00087                                         GiNaC::ex Vn = inner(V, normal_vec);
00088                                         Vn = Vn.subs(line.repr(s).op(0)).subs(line.repr(s).op(1));
00089                                         basis2coeff = pol2basisandcoeff(Vn,s);
00090                                         for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
00091                                         {
00092                                                 GiNaC::ex basis = (*iter).first;
00093                                                 GiNaC::ex coeff= (*iter).second;
00094                                                 if ( coeff != 0 &&  basis.degree(s) > int(order+1) )
00095                                                 {
00096                                                         equations.append( coeff == 0 );
00097                                                 }
00098                                         }
00099                                 }
00100 
00101                                 if (  order%2==1  )
00102                                 {
00103                                         // tangent constraints on edges:
00104                                         for (int i=0; i< 3; i++)
00105                                         {
00106                                                 Line line = triangle.line(i);
00107                                                 GiNaC::symbol s("s");
00108                                                 GiNaC::lst tangent_vec = tangent(triangle, i);
00109                                                 GiNaC::ex Vt = inner(V, tangent_vec);
00110                                                 Vt = Vt.subs(line.repr(s).op(0)).subs(line.repr(s).op(1));
00111                                                 basis2coeff = pol2basisandcoeff(Vt,s);
00112                                                 for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
00113                                                 {
00114                                                         GiNaC::ex basis = (*iter).first;
00115                                                         GiNaC::ex coeff= (*iter).second;
00116                                                         if ( coeff != 0 &&  basis.degree(s) > int(order+2) )
00117                                                         {
00118                                                                 equations.append( coeff == 0 );
00119                                                         }
00120                                                 }
00121                                         }
00122                                 }
00123 
00124                                 GiNaC::ex dofi;
00125 
00126                                 // dofs related to the normal on the edges
00127                                 for (int i=0; i< 3; i++)
00128                                 {
00129                                         Line line = triangle.line(i);
00130                                         GiNaC::lst normal_vec = normal(triangle, i);
00131                                         GiNaC::ex Vn = inner(V, normal_vec);
00132                                         GiNaC::lst points = interior_coordinates(line, order + 1);
00133                                         GiNaC::ex edge_length = line.integrate(GiNaC::numeric(1));
00134 
00135                                         GiNaC::ex point;
00136                                         for (unsigned int j=0; j< points.nops(); j++)
00137                                         {
00138                                                 point = points.op(j);
00139                                                 dofi = Vn.subs(x == point.op(0)).subs(y == point.op(1))*edge_length;
00140                                                 dofs.insert(dofs.end(), GiNaC::lst(point,nsymb));
00141                                                 GiNaC::ex eq = dofi == GiNaC::numeric(0);
00142                                                 equations.append(eq);
00143                                         }
00144                                 }
00145 
00146                                 if ( order%2==0)
00147                                 {
00148                                         // dofs related to the tangent on the edges
00149                                         for (int i=0; i< 3; i++)
00150                                         {
00151                                                 Line line = triangle.line(i);
00152                                                 GiNaC::lst tangent_vec = tangent(triangle, i);
00153                                                 GiNaC::ex Vt = inner(V, tangent_vec);
00154                                                 GiNaC::lst points = interior_coordinates(line, order);
00155                                                 GiNaC::ex edge_length = line.integrate(GiNaC::numeric(1));
00156 
00157                                                 GiNaC::ex point;
00158                                                 for (unsigned int j=0; j< points.nops(); j++)
00159                                                 {
00160                                                         point = points.op(j);
00161                                                         dofi = Vt.subs(x == point.op(0)).subs(y == point.op(1))*edge_length;
00162                                                         dofs.insert(dofs.end(), GiNaC::lst(point,tsymb));
00163                                                         GiNaC::ex eq = dofi == GiNaC::numeric(0);
00164                                                         equations.append(eq);
00165                                                 }
00166                                         }
00167                                 }
00168 
00169                                 if ( order%2==1 )
00170                                 {
00171                                         // dofs related to the tangent on the edges
00172                                         for (int i=0; i< 3; i++)
00173                                         {
00174                                                 Line line = triangle.line(i);
00175                                                 GiNaC::lst tangent_vec = tangent(triangle, i);
00176                                                 GiNaC::ex Vt = inner(V, tangent_vec);
00177                                                 GiNaC::lst points = interior_coordinates(line, order-1);
00178                                                 GiNaC::ex edge_length = line.integrate(GiNaC::numeric(1));
00179 
00180                                                 GiNaC::ex point;
00181                                                 for (unsigned int j=0; j< points.nops(); j++)
00182                                                 {
00183                                                         point = points.op(j);
00184                                                         dofi = Vt.subs(x == point.op(0)).subs(y == point.op(1))*edge_length;
00185                                                         dofs.insert(dofs.end(), GiNaC::lst(point,tsymb));
00186                                                         GiNaC::ex eq = dofi == GiNaC::numeric(0);
00187                                                         equations.append(eq);
00188                                                 }
00189                                         }
00190                                 }
00191 
00192                                 // dofs related to the whole triangle
00193                                 GiNaC::lst bernstein_polv;
00194                                 if ( order > 0)
00195                                 {
00196                                         GiNaC::lst points = interior_coordinates(triangle, order-1);
00197                                         GiNaC::ex point;
00198                                         GiNaC::ex eq;
00199                                         for (unsigned int i=0; i< points.nops(); i++)
00200                                         {
00201                                                 point = points.op(i);
00202 
00203                                                 // x - components of interior dofs
00204                                                 dofi = V.op(0).subs(x == point.op(0)).subs(y == point.op(1));
00205                                                 eq = dofi == GiNaC::numeric(0);
00206                                                 equations.append(eq);
00207                                                 dofs.insert(dofs.end(), GiNaC::lst(point,0));
00208 
00209                                                 // y - components of interior dofs
00210                                                 dofi = V.op(1).subs(x == point.op(0)).subs(y == point.op(1));
00211                                                 eq = dofi == GiNaC::numeric(0);
00212                                                 equations.append(eq);
00213                                                 dofs.insert(dofs.end(), GiNaC::lst(point,1));
00214                                         }
00215                                 }
00216 
00217                                 variables = collapse(GiNaC::lst(V_vars.op(0),V_vars.op(1)));
00218 
00219                                 //              std::cout <<"no variables "<<variables.nops()<<std::endl;
00220                                 //              std::cout <<"no equations "<<equations.nops()<<std::endl;
00221 
00222                                 GiNaC::matrix b; GiNaC::matrix A;
00223                                 matrix_from_equations(equations, variables, A, b);
00224 
00225                                 //      std::cout <<" A "<<A.evalf()<<std::endl;
00226 
00227                                 unsigned int ncols = A.cols();
00228                                 GiNaC::matrix vars_sq(ncols, ncols);
00229 
00230                                 // matrix of symbols
00231                                 for (unsigned r=0; r<ncols; ++r)
00232                                         for (unsigned c=0; c<ncols; ++c)
00233                                                 vars_sq(r, c) = GiNaC::symbol();
00234 
00235                                 GiNaC::matrix id(ncols, ncols);
00236 
00237                                 // identity
00238                                 const GiNaC::ex _ex1(1);
00239                                 for (unsigned i=0; i<ncols; ++i)
00240                                         id(i, i) = _ex1;
00241 
00242                                 // invert the matrix
00243                                 GiNaC::matrix m_inv(ncols, ncols);
00244                                 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss);
00245 
00246                                 for (unsigned int i=0; i<dofs.size(); i++)
00247                                 {
00248                                         b.let_op(11+i) = GiNaC::numeric(1);
00249                                         GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b));
00250 
00251                                         GiNaC::lst subs;
00252                                         for (unsigned int ii=0; ii<xx.nops(); ii++)
00253                                         {
00254                                                 subs.append(variables.op(ii) == xx.op(ii));
00255                                         }
00256                                         GiNaC::ex Nj1 = V.op(0).subs(subs);
00257                                         GiNaC::ex Nj2 = V.op(1).subs(subs);
00258                                         Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2)));
00259                                         b.let_op(11+i) = GiNaC::numeric(0);
00260 
00261                                 }
00262                         }                                        //dofs are integrals etc. ie not represented as normal/tangents in points
00263                         else
00264                         {
00265 
00266                                 description = "Robust_2D";
00267 
00268                                 Triangle& triangle = (Triangle&)(*p);
00269                                 GiNaC::lst equations;
00270                                 GiNaC::lst variables;
00271                                 GiNaC::ex V_space = bernsteinv(2, order+3, triangle, "a");
00272                                 GiNaC::ex V = V_space.op(0);
00273                                 GiNaC::ex V_vars = V_space.op(1);
00274 
00275                                 GiNaC::ex divV = div(V);
00276                                 exmap basis2coeff = pol2basisandcoeff(divV);
00277                                 exmap::iterator iter;
00278 
00279                                 // div constraints:
00280                                 for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
00281                                 {
00282                                         GiNaC::ex basis = (*iter).first;
00283                                         GiNaC::ex coeff= (*iter).second;
00284                                         if ( coeff != 0 )
00285                                         {
00286                                         }
00287                                         if ( coeff != 0 && ( basis.degree(x) + basis.degree(y)  > int(order) )  )
00288                                         {
00289                                                 equations.append( coeff == 0 );
00290                                         }
00291                                 }
00292 
00293                                 // normal constraints on edges:
00294                                 for (int i=0; i< 3; i++)
00295                                 {
00296                                         Line line = triangle.line(i);
00297                                         GiNaC::symbol s("s");
00298                                         GiNaC::lst normal_vec = normal(triangle, i);
00299                                         GiNaC::ex Vn = inner(V, normal_vec);
00300                                         Vn = Vn.subs(line.repr(s).op(0)).subs(line.repr(s).op(1));
00301                                         basis2coeff = pol2basisandcoeff(Vn,s);
00302                                         for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
00303                                         {
00304                                                 GiNaC::ex basis = (*iter).first;
00305                                                 GiNaC::ex coeff= (*iter).second;
00306                                                 if ( coeff != 0 &&  basis.degree(s) > int(order+1) )
00307                                                 {
00308                                                         equations.append( coeff == 0 );
00309                                                 }
00310                                         }
00311                                 }
00312 
00313                                 if (  order%2==1  )
00314                                 {
00315                                         // tangent constraints on edges:
00316                                         for (int i=0; i< 3; i++)
00317                                         {
00318                                                 Line line = triangle.line(i);
00319                                                 GiNaC::symbol s("s");
00320                                                 GiNaC::lst tangent_vec = tangent(triangle, i);
00321                                                 GiNaC::ex Vt = inner(V, tangent_vec);
00322                                                 Vt = Vt.subs(line.repr(s).op(0)).subs(line.repr(s).op(1));
00323                                                 basis2coeff = pol2basisandcoeff(Vt,s);
00324                                                 for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
00325                                                 {
00326                                                         GiNaC::ex basis = (*iter).first;
00327                                                         GiNaC::ex coeff= (*iter).second;
00328                                                         if ( coeff != 0 &&  basis.degree(s) > int(order+2) )
00329                                                         {
00330                                                                 equations.append( coeff == 0 );
00331                                                         }
00332                                                 }
00333                                         }
00334                                 }
00335 
00336                                 // dofs related to the normal on the edges
00337                                 for (int i=0; i< 3; i++)
00338                                 {
00339                                         Line line = triangle.line(i);
00340                                         GiNaC::lst normal_vec = normal(triangle, i);
00341                                         GiNaC::ex Pk1_space = bernstein(order+1, line, istr("a",i));
00342                                         GiNaC::ex Pk1 = Pk1_space.op(2);
00343                                         GiNaC::ex Vn = inner(V, normal_vec);
00344 
00345                                         GiNaC::ex basis;
00346                                         for (unsigned int j=0; j< Pk1.nops(); j++)
00347                                         {
00348                                                 basis = Pk1.op(j);
00349                                                 GiNaC::ex integrand = Vn*basis;
00350                                                 GiNaC::ex dofi =  line.integrate(integrand);
00351                                                 //                dofs.insert(dofs.end(), dofi);
00352                                                 GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),j);
00353                                                 dofs.insert(dofs.end(), d);
00354                                                 GiNaC::ex eq = dofi == GiNaC::numeric(0);
00355                                                 equations.append(eq);
00356 
00357                                                 GiNaC::ex u = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]")));
00358                                                 GiNaC::ex n = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("normal_vec[0]"), GiNaC::symbol("normal_vec[1]")));
00359                                                 dof_repr.append(GiNaC::lst(inner(u,n)*basis.subs( x == GiNaC::symbol("xi[0]"))
00360                                                         .subs( y == GiNaC::symbol("xi[1]")), d));
00361 
00362                                         }
00363                                 }
00364 
00365                                 // dofs related to the tangent on the edges
00366                                 if (  order%2==0  )
00367                                 {
00368                                         for (int i=0; i< 3; i++)
00369                                         {
00370                                                 Line line = triangle.line(i);
00371                                                 GiNaC::lst tangent_vec = tangent(triangle, i);
00372                                                 GiNaC::ex Pk_space = bernstein(order, line, istr("a",i));
00373                                                 GiNaC::ex Pk = Pk_space.op(2);
00374                                                 GiNaC::ex Vt = inner(V, tangent_vec);
00375 
00376                                                 GiNaC::ex basis;
00377                                                 for (unsigned int j=0; j< Pk.nops(); j++)
00378                                                 {
00379                                                         basis = Pk.op(j);
00380                                                         GiNaC::ex integrand = Vt*basis;
00381                                                         GiNaC::ex dofi =  line.integrate(integrand);
00382                                                         //                dofs.insert(dofs.end(), dofi);
00383                                                         GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),2);
00384                                                         dofs.insert(dofs.end(), d);
00385                                                         GiNaC::ex eq = dofi == GiNaC::numeric(0);
00386                                                         equations.append(eq);
00387 
00388                                                         GiNaC::ex u = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]")));
00389                                                         GiNaC::ex t = GiNaC::matrix(2,1,GiNaC::lst(-GiNaC::symbol("normal_vec[1]"), GiNaC::symbol("normal_vec[0]")));
00390                                                         dof_repr.append(GiNaC::lst(inner(u,t)*basis.subs( x == GiNaC::symbol("xi[0]"))
00391                                                                 .subs( y == GiNaC::symbol("xi[1]")), d));
00392 
00393                                                 }
00394                                         }
00395                                 }
00396 
00397                                 if (  order%2==1  )
00398                                 {
00399                                         for (int i=0; i< 3; i++)
00400                                         {
00401                                                 Line line = triangle.line(i);
00402                                                 GiNaC::lst tangent_vec = tangent(triangle, i);
00403                                                 GiNaC::ex Pk_space = bernstein(order-1, line, istr("a",i));
00404                                                 GiNaC::ex Pk = Pk_space.op(2);
00405                                                 GiNaC::ex Vt = inner(V, tangent_vec);
00406 
00407                                                 GiNaC::ex basis;
00408                                                 for (unsigned int j=0; j< Pk.nops(); j++)
00409                                                 {
00410                                                         basis = Pk.op(j);
00411                                                         GiNaC::ex integrand = Vt*basis;
00412                                                         GiNaC::ex dofi =  line.integrate(integrand);
00413                                                         //                dofs.insert(dofs.end(), dofi);
00414                                                         GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),2);
00415                                                         dofs.insert(dofs.end(), d);
00416                                                         GiNaC::ex eq = dofi == GiNaC::numeric(0);
00417                                                         equations.append(eq);
00418 
00419                                                         GiNaC::ex u = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]")));
00420                                                         GiNaC::ex t = GiNaC::matrix(2,1,GiNaC::lst(-GiNaC::symbol("normal_vec[1]"), GiNaC::symbol("normal_vec[0]")));
00421                                                         dof_repr.append(GiNaC::lst(inner(u,t)*basis.subs( x == GiNaC::symbol("xi[0]"))
00422                                                                 .subs( y == GiNaC::symbol("xi[1]")), d));
00423 
00424                                                 }
00425                                         }
00426                                 }
00427 
00428                                 // dofs related to the whole triangle
00429                                 GiNaC::lst bernstein_polv;
00430                                 if ( order > 0)
00431                                 {
00432                                         bernstein_polv = bernsteinv(2,order-1, triangle, "a");
00433                                         GiNaC::ex basis_space = bernstein_polv.op(2);
00434                                         for (unsigned int i=0; i< basis_space.nops(); i++)
00435                                         {
00436                                                 GiNaC::lst basis = GiNaC::ex_to<GiNaC::lst>(basis_space.op(i));
00437                                                 GiNaC::ex integrand = inner(V, basis);
00438                                                 GiNaC::ex dofi = triangle.integrate(integrand);
00439                                                 GiNaC::ex eq = dofi == GiNaC::numeric(0);
00440                                                 equations.append(eq);
00441 
00442                                                 GiNaC::lst d = GiNaC::lst(GiNaC::lst(triangle.vertex(0), triangle.vertex(1), triangle.vertex(2)), i);
00443                                                 dofs.insert(dofs.end(), d);
00444 
00445                                         }
00446                                 }
00447 
00448                                 variables = collapse(GiNaC::lst(V_vars.op(0),V_vars.op(1)));
00449 
00450                                 //              std::cout <<"no variables "<<variables.nops()<<std::endl;
00451                                 //              std::cout <<"no equations "<<equations.nops()<<std::endl;
00452 
00453                                 GiNaC::matrix b; GiNaC::matrix A;
00454                                 matrix_from_equations(equations, variables, A, b);
00455 
00456                                 //        std::cout <<" A "<<A.evalf()<<std::endl;
00457 
00458                                 unsigned int ncols = A.cols();
00459                                 GiNaC::matrix vars_sq(ncols, ncols);
00460 
00461                                 // matrix of symbols
00462                                 for (unsigned r=0; r<ncols; ++r)
00463                                         for (unsigned c=0; c<ncols; ++c)
00464                                                 vars_sq(r, c) = GiNaC::symbol();
00465 
00466                                 GiNaC::matrix id(ncols, ncols);
00467 
00468                                 // identity
00469                                 const GiNaC::ex _ex1(1);
00470                                 for (unsigned i=0; i<ncols; ++i)
00471                                         id(i, i) = _ex1;
00472 
00473                                 // invert the matrix
00474                                 GiNaC::matrix m_inv(ncols, ncols);
00475                                 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss);
00476 
00477                                 for (unsigned int i=0; i<dofs.size(); i++)
00478                                 {
00479                                         b.let_op(11+i) = GiNaC::numeric(1);
00480                                         GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b));
00481 
00482                                         GiNaC::lst subs;
00483                                         for (unsigned int ii=0; ii<xx.nops(); ii++)
00484                                         {
00485                                                 subs.append(variables.op(ii) == xx.op(ii));
00486                                         }
00487                                         GiNaC::ex Nj1 = V.op(0).subs(subs);
00488                                         GiNaC::ex Nj2 = V.op(1).subs(subs);
00489                                         Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2)));
00490                                         b.let_op(11+i) = GiNaC::numeric(0);
00491                                 }
00492                         }
00493                 }
00494         }
00495 
00496         void Robust:: compute_basis_functions_old()
00497         {
00498 
00499                 // remove previously computed basis functions and dofs
00500                 Ns.clear();
00501                 dofs.clear();
00502 
00503                 order = 3;
00504 
00505                 if ( p == NULL )
00506                 {
00507                         throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00508                 }
00509 
00510                 if ( p->str().find("Line") != string::npos )
00511                 {
00512 
00513                         cout <<"Can not define the Robust element on a line"<<endl;
00514 
00515                 }
00516                 else if ( p->str().find("Triangle") != string::npos )
00517                 {
00518 
00519                         description = "Robust_2D";
00520 
00521                         Triangle& triangle = (Triangle&)(*p);
00522                         GiNaC::lst equations;
00523                         GiNaC::lst variables;
00524                         GiNaC::ex V_space = bernsteinv(2, order, triangle, "a");
00525                         GiNaC::ex V = V_space.op(0);
00526                         GiNaC::ex V_vars = V_space.op(1);
00527 
00528                         GiNaC::ex divV = div(V);
00529                         exmap basis2coeff = pol2basisandcoeff(divV);
00530                         exmap::iterator iter;
00531 
00532                         // div constraints:
00533                         for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
00534                         {
00535                                 GiNaC::ex basis = (*iter).first;
00536                                 GiNaC::ex coeff= (*iter).second;
00537                                 if ( coeff != 0 && ( basis.degree(x) > 0 || basis.degree(y)  > 0 )  )
00538                                 {
00539                                         equations.append( coeff == 0 );
00540                                 }
00541                         }
00542 
00543                         // constraints on edges:
00544                         for (int i=0; i< 3; i++)
00545                         {
00546                                 Line line = triangle.line(i);
00547                                 GiNaC::symbol s("s");
00548                                 GiNaC::lst normal_vec = normal(triangle, i);
00549                                 GiNaC::ex Vn = inner(V, normal_vec);
00550                                 Vn = Vn.subs(line.repr(s).op(0)).subs(line.repr(s).op(1));
00551                                 basis2coeff = pol2basisandcoeff(Vn,s);
00552                                 for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
00553                                 {
00554                                         GiNaC::ex basis = (*iter).first;
00555                                         GiNaC::ex coeff= (*iter).second;
00556                                         if ( coeff != 0 &&  basis.degree(s) > 1 )
00557                                         {
00558                                                 equations.append( coeff == 0 );
00559                                         }
00560                                 }
00561                         }
00562 
00563                         // dofs related to the normal on the edges
00564                         for (int i=0; i< 3; i++)
00565                         {
00566                                 Line line = triangle.line(i);
00567                                 GiNaC::lst normal_vec = normal(triangle, i);
00568                                 GiNaC::ex Pk1_space = bernstein(1, line, istr("a",i));
00569                                 GiNaC::ex Pk1 = Pk1_space.op(2);
00570                                 GiNaC::ex Vn = inner(V, normal_vec);
00571 
00572                                 GiNaC::ex basis;
00573                                 for (unsigned int j=0; j< Pk1.nops(); j++)
00574                                 {
00575                                         basis = Pk1.op(j);
00576                                         GiNaC::ex integrand = Vn*basis;
00577                                         GiNaC::ex dofi =  line.integrate(integrand);
00578                                         //                dofs.insert(dofs.end(), dofi);
00579                                         GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),j);
00580                                         dofs.insert(dofs.end(), d);
00581                                         GiNaC::ex eq = dofi == GiNaC::numeric(0);
00582                                         equations.append(eq);
00583 
00584                                         GiNaC::ex u = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]")));
00585                                         GiNaC::ex n = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("normal_vec[0]"), GiNaC::symbol("normal_vec[1]")));
00586                                         dof_repr.append(GiNaC::lst(inner(u,n)*basis.subs( x == GiNaC::symbol("xi[0]"))
00587                                                 .subs( y == GiNaC::symbol("xi[1]")), d));
00588 
00589                                 }
00590                         }
00591 
00592                         // dofs related to the tangent on the edges
00593                         for (int i=0; i< 3; i++)
00594                         {
00595                                 Line line = triangle.line(i);
00596                                 GiNaC::lst tangent_vec = tangent(triangle, i);
00597                                 GiNaC::ex Pk_space = bernstein(0, line, istr("a",i));
00598                                 GiNaC::ex Pk = Pk_space.op(2);
00599                                 GiNaC::ex Vt = inner(V, tangent_vec);
00600 
00601                                 GiNaC::ex basis;
00602                                 for (unsigned int j=0; j< Pk.nops(); j++)
00603                                 {
00604                                         basis = Pk.op(j);
00605                                         GiNaC::ex integrand = Vt*basis;
00606                                         GiNaC::ex dofi =  line.integrate(integrand);
00607                                         //                dofs.insert(dofs.end(), dofi);
00608                                         GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),2);
00609                                         dofs.insert(dofs.end(), d);
00610                                         GiNaC::ex eq = dofi == GiNaC::numeric(0);
00611                                         equations.append(eq);
00612 
00613                                         GiNaC::ex u = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]")));
00614                                         GiNaC::ex t = GiNaC::matrix(2,1,GiNaC::lst(-GiNaC::symbol("normal_vec[1]"), GiNaC::symbol("normal_vec[0]")));
00615                                         dof_repr.append(GiNaC::lst(inner(u,t)*basis.subs( x == GiNaC::symbol("xi[0]"))
00616                                                 .subs( y == GiNaC::symbol("xi[1]")), d));
00617 
00618                                 }
00619                         }
00620 
00621                         variables = collapse(GiNaC::lst(V_vars.op(0),V_vars.op(1)));
00622 
00623                         GiNaC::matrix b; GiNaC::matrix A;
00624                         matrix_from_equations(equations, variables, A, b);
00625 
00626                         unsigned int ncols = A.cols();
00627                         GiNaC::matrix vars_sq(ncols, ncols);
00628 
00629                         // matrix of symbols
00630                         for (unsigned r=0; r<ncols; ++r)
00631                                 for (unsigned c=0; c<ncols; ++c)
00632                                         vars_sq(r, c) = GiNaC::symbol();
00633 
00634                         GiNaC::matrix id(ncols, ncols);
00635 
00636                         // identity
00637                         const GiNaC::ex _ex1(1);
00638                         for (unsigned i=0; i<ncols; ++i)
00639                                 id(i, i) = _ex1;
00640 
00641                         // invert the matrix
00642                         GiNaC::matrix m_inv(ncols, ncols);
00643                         m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss);
00644 
00645                         for (unsigned int i=0; i<dofs.size(); i++)
00646                         {
00647                                 b.let_op(11+i) = GiNaC::numeric(1);
00648                                 GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b));
00649 
00650                                 GiNaC::lst subs;
00651                                 for (unsigned int ii=0; ii<xx.nops(); ii++)
00652                                 {
00653                                         subs.append(variables.op(ii) == xx.op(ii));
00654                                 }
00655                                 GiNaC::ex Nj1 = V.op(0).subs(subs);
00656                                 GiNaC::ex Nj2 = V.op(1).subs(subs);
00657                                 Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2)));
00658                                 b.let_op(11+i) = GiNaC::numeric(0);
00659                         }
00660                 }
00661         }
00662 
00663 }                                                                // namespace SyFi
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines