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