SyFi 0.3
|
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