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

Generated on Mon Aug 31 16:16:45 2009 for SyFi by  doxygen 1.5.9