00001
00002
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
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
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
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
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
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
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
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
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
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
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
00220
00221
00222 GiNaC::matrix b; GiNaC::matrix A;
00223 matrix_from_equations(equations, variables, A, b);
00224
00225
00226
00227 unsigned int ncols = A.cols();
00228 GiNaC::matrix vars_sq(ncols, ncols);
00229
00230
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
00238 const GiNaC::ex _ex1(1);
00239 for (unsigned i=0; i<ncols; ++i)
00240 id(i, i) = _ex1;
00241
00242
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 }
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
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
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
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
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
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
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
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
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
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
00451
00452
00453 GiNaC::matrix b; GiNaC::matrix A;
00454 matrix_from_equations(equations, variables, A, b);
00455
00456
00457
00458 unsigned int ncols = A.cols();
00459 GiNaC::matrix vars_sq(ncols, ncols);
00460
00461
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
00469 const GiNaC::ex _ex1(1);
00470 for (unsigned i=0; i<ncols; ++i)
00471 id(i, i) = _ex1;
00472
00473
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
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
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
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
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
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
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
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
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
00637 const GiNaC::ex _ex1(1);
00638 for (unsigned i=0; i<ncols; ++i)
00639 id(i, i) = _ex1;
00640
00641
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 }