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