00001
00002
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
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
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
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
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
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
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
00153
00154
00155
00156
00157
00158
00159
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
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
00175 const GiNaC::ex _ex1(1);
00176 for (unsigned i=0; i<ncols; ++i)
00177 id(i, i) = _ex1;
00178
00179
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
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
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
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
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
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
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
00306
00307
00308
00309
00310
00311
00312
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
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
00328 const GiNaC::ex _ex1(1);
00329 for (unsigned i=0; i<ncols; ++i)
00330 id(i, i) = _ex1;
00331
00332
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
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
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
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
00467
00468
00469
00470
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
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
00486 const GiNaC::ex _ex1(1);
00487 for (unsigned i=0; i<ncols; ++i)
00488 id(i, i) = _ex1;
00489
00490
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
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
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
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
00618
00619
00620
00621
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(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 }