|
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 <sstream> 00019 00020 #include "Lagrange.h" 00021 #include "ElementComputations.h" 00022 #include "tools.h" 00023 00024 using std::cout; 00025 using std::endl; 00026 using std::string; 00027 00028 namespace SyFi 00029 { 00030 00031 Lagrange::Lagrange() : StandardFE() 00032 { 00033 description = "Lagrange"; 00034 } 00035 00036 Lagrange::Lagrange(Polygon& p, unsigned int order) : StandardFE(p, order) 00037 { 00038 compute_basis_functions(); 00039 } 00040 00041 void Lagrange:: compute_basis_functions() 00042 { 00043 00044 // NOTE: in the below code dof(i) is not used to 00045 // determine the basis functions 00046 00047 // remove previously computed basis functions and dofs 00048 Ns.clear(); 00049 dofs.clear(); 00050 00051 if ( order < 1 ) 00052 { 00053 throw(std::logic_error("Lagrangian elements must be of order 1 or higher.")); 00054 } 00055 00056 if ( p == NULL ) 00057 { 00058 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed")); 00059 } 00060 00061 GiNaC::lst equations; 00062 GiNaC::lst variables; 00063 GiNaC::ex polynom; 00064 00065 if (p->str().find("ReferenceLine") != string::npos) 00066 { 00067 00068 description = istr("Lagrange_", order) + "_1D"; 00069 00070 // Look at the case with the Triangle for a documented code 00071 00072 // polynom = pol(order, 1, "a"); 00073 // variables = coeffs(polynom); 00074 GiNaC::ex polynom_space = bernstein(order, *p, "a"); 00075 // GiNaC::ex polynom_space = pol(order, 1, "a"); 00076 polynom = polynom_space.op(0); 00077 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00078 00079 GiNaC::ex increment = GiNaC::numeric(1,order); 00080 00081 GiNaC::ex Nj; 00082 for (GiNaC::ex p=0; p<= 1 ; p += increment ) 00083 { 00084 GiNaC::ex eq = polynom == GiNaC::numeric(0); 00085 equations.append(eq.subs(x == p)); 00086 dofs.insert(dofs.end(), GiNaC::lst(p)); 00087 } 00088 } 00089 else if (p->str().find("Line") != string::npos ) 00090 { 00091 description = istr("Lagrange_", order) + "_1D"; 00092 00093 // Look at the case with the Triangle for a documented code 00094 00095 // polynom = pol(order, 1, "a"); 00096 // variables = coeffs(polynom); 00097 GiNaC::ex polynom_space = bernstein(order, *p, "a"); 00098 // GiNaC::ex polynom_space = pol(order, 1, "a"); 00099 polynom = polynom_space.op(0); 00100 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00101 00102 Polygon& pp = *p; 00103 Line& l = (Line&) pp; 00104 GiNaC::lst points = bezier_ordinates(l,order); 00105 00106 GiNaC::ex Nj; 00107 for (unsigned int i=1; i<= points.nops() ; i++ ) 00108 { 00109 GiNaC::ex point = points.op(i-1); 00110 GiNaC::ex eq = polynom == GiNaC::numeric(0); 00111 if (point.nops() == 0) eq = eq.subs(x == point); 00112 if (point.nops() > 0) eq = eq.subs(x == point.op(0)); 00113 if (point.nops() > 1) eq = eq.subs(y == point.op(1)); 00114 if (point.nops() > 2) eq = eq.subs(z == point.op(2)); 00115 equations.append(eq); 00116 dofs.insert(dofs.end(), GiNaC::lst(point)); 00117 } 00118 00119 } 00120 else if (p->str().find("ReferenceTriangle") != string::npos ) 00121 { 00122 00123 description = istr("Lagrange_", order) + "_2D"; 00124 00125 // Look at the case with the Triangle for a documented code 00126 00127 // polynom = pol(order, 2, "b"); 00128 // variables = coeffs(polynom); 00129 GiNaC::ex polynom_space = bernstein(order, *p, "b"); 00130 // GiNaC::ex polynom_space = pol(order, 2, "a"); 00131 polynom = polynom_space.op(0); 00132 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00133 00134 GiNaC::ex increment = GiNaC::numeric(1,order); 00135 00136 GiNaC::ex Nj; 00137 GiNaC::numeric one = 1; 00138 for (GiNaC::ex q=0; q<= one ; q += increment ) 00139 { 00140 for (GiNaC::ex p=0; p<= one-q ; p += increment ) 00141 { 00142 GiNaC::ex eq = polynom == GiNaC::numeric(0); 00143 equations.append(eq.subs(GiNaC::lst(x == p, y == q))); 00144 dofs.insert(dofs.end(), GiNaC::lst(p,q)); 00145 } 00146 } 00147 } 00148 else if ( p->str().find("Triangle") != string::npos) 00149 { 00150 00151 description = istr("Lagrange_", order) + "_2D"; 00152 00153 // Look HERE for the documented code 00154 00155 // GiNaC::ex polynom_space = pol(order, 2, "a"); 00156 GiNaC::ex polynom_space = bernstein(order, *p, "b"); 00157 // the polynomial spaces on the form: 00158 // first item: a0 + a1*x + a2*y + a3*x^2 + a4*x*y ... the polynom 00159 // second item: a0, a1, a2, ... the coefficents 00160 // third item 1, x, y, x^2, .. the basis 00161 // Could also do: 00162 // GiNaC::ex polynom_space = bernstein(order, t, "a"); 00163 00164 polynom = polynom_space.op(0); 00165 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00166 00167 GiNaC::ex Nj; 00168 Polygon& pp = *p; 00169 Triangle& t = (Triangle&) pp; 00170 // The bezier ordinates (in which the basis function should be either 0 or 1) 00171 GiNaC::lst points = bezier_ordinates(t,order); 00172 00173 for (unsigned int i=1; i<= points.nops() ; i++ ) 00174 { 00175 GiNaC::ex point = points.op(i-1); 00176 GiNaC::ex eq = polynom == GiNaC::numeric(0); 00177 equations.append(eq.subs(GiNaC::lst(x == point.op(0) , y == point.op(1)))); 00178 dofs.insert(dofs.end(), GiNaC::lst(point.op(0),point.op(1))); 00179 } 00180 } 00181 00182 else if ( p->str().find("ReferenceRectangle") != string::npos) 00183 { 00184 00185 description = istr("Lagrange_", order) + "_2D"; 00186 00187 // create 1D element, then create tensor product 00188 ReferenceLine line; 00189 Lagrange fe(line, order); 00190 00191 for (unsigned int i=0; i< fe.nbf(); i++) 00192 { 00193 for (unsigned int j=0; j< fe.nbf(); j++) 00194 { 00195 Ns.insert(Ns.end(), fe.N(i)*fe.N(j).subs(x==y)); 00196 dofs.insert(dofs.end(), GiNaC::lst(fe.dof(i).op(0), fe.dof(j).op(0))); 00197 } 00198 } 00199 return; 00200 00201 /* OLD CODE 00202 00203 GiNaC::ex polynom_space = legendre(order, 2, "b"); 00204 00205 polynom = polynom_space.op(0); 00206 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00207 00208 GiNaC::ex Nj; 00209 Polygon& pp = *p; 00210 ReferenceRectangle& b = (ReferenceRectangle&) pp; 00211 00212 int no_points = (order+1)*(order+1); 00213 GiNaC::ex increment = GiNaC::numeric(1,order); 00214 00215 GiNaC::numeric one=1.0; 00216 00217 for (GiNaC::ex q=0; q <= one ; q += increment ) 00218 { 00219 for (GiNaC::ex p=0; p <= one ; p += increment ) 00220 { 00221 GiNaC::ex eq = polynom == GiNaC::numeric(0); 00222 equations.append(eq.subs(GiNaC::lst(x == p, y == q))); 00223 dofs.push_back(GiNaC::lst(p,q)); 00224 } 00225 } 00226 */ 00227 } 00228 else if ( p->str().find("ReferenceTetrahedron") != string::npos) 00229 { 00230 00231 description = istr("Lagrange_", order) + "_3D"; 00232 00233 // Look at the case with the Triangle for a documented code 00234 00235 // polynom = pol(order, 3, "b"); 00236 // GiNaC::ex polynom_space = pol(order, 3, "a"); 00237 GiNaC::ex polynom_space = bernstein(order, *p, "b"); 00238 polynom = polynom_space.op(0); 00239 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00240 00241 int nno =0; 00242 for (unsigned int j=0; j<= order; j++) 00243 { 00244 nno += (j+1)*(j+2)/2; 00245 } 00246 00247 GiNaC::ex increment = GiNaC::numeric(1,order); 00248 00249 GiNaC::ex Nj; 00250 for (GiNaC::ex r=0; r<= 1 ; r += increment ) 00251 { 00252 for (GiNaC::ex q=0; q<= 1-r ; q += increment ) 00253 { 00254 for (GiNaC::ex p=0; p<= 1-r-q ; p += increment ) 00255 { 00256 GiNaC::ex eq = polynom == GiNaC::numeric(0); 00257 equations.append(eq.subs(GiNaC::lst(x == p, y == q, z == r ))); 00258 dofs.push_back(GiNaC::lst(p,q,r)); 00259 } 00260 } 00261 } 00262 } 00263 else if ( p->str().find("Tetrahedron") != string::npos) 00264 { 00265 00266 description = istr("Lagrange_", order) + "_3D"; 00267 00268 // Look at the case with the Triangle for a documented code 00269 00270 GiNaC::ex polynom_space = pol(order, 3, "a"); 00271 // GiNaC::ex polynom_space = bernstein(order, *p, "b"); 00272 polynom = polynom_space.op(0); 00273 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00274 00275 GiNaC::ex increment = GiNaC::numeric(1,order); 00276 00277 GiNaC::ex Nj; 00278 Polygon& pp = *p; 00279 Tetrahedron& t = (Tetrahedron&) pp; 00280 GiNaC::lst points = bezier_ordinates(t,order); 00281 for (unsigned int i=1; i<= points.nops() ; i++ ) 00282 { 00283 GiNaC::ex point = points.op(i-1); 00284 GiNaC::ex eq = polynom == GiNaC::numeric(0); 00285 equations.append(eq.subs(GiNaC::lst(x == point.op(0) , y == point.op(1), z == point.op(2)))); 00286 dofs.push_back(GiNaC::lst(point.op(0),point.op(1),point.op(2))); 00287 } 00288 } 00289 else if ( p->str().find("ReferenceBox") != string::npos) 00290 { 00291 00292 description = istr("Lagrange_", order) + "_3D"; 00293 00294 ReferenceLine line; 00295 Lagrange fe(line, order); 00296 00297 for (unsigned int i=0; i< fe.nbf(); i++) 00298 { 00299 for (unsigned int j=0; j< fe.nbf(); j++) 00300 { 00301 for (unsigned int k=0; k< fe.nbf(); k++) 00302 { 00303 Ns.insert(Ns.end(), fe.N(i)*fe.N(j).subs(x==y)*fe.N(k).subs(x==z)); 00304 dofs.insert(dofs.end(), GiNaC::lst(fe.dof(i).op(0), fe.dof(j).op(0), fe.dof(k).op(0))); 00305 } 00306 } 00307 } 00308 return; 00309 00310 /* OLD CODE 00311 00312 GiNaC::ex polynom_space = legendre(order, 3, "b"); 00313 00314 polynom = polynom_space.op(0); 00315 variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)); 00316 00317 GiNaC::ex Nj; 00318 Polygon& pp = *p; 00319 ReferenceRectangle& b = (ReferenceRectangle&) pp; 00320 00321 int no_points = (order+1)*(order+1)*(order+1); 00322 GiNaC::ex increment = GiNaC::numeric(1,order); 00323 00324 GiNaC::numeric one = 1; 00325 for (GiNaC::ex p=0; p <= one ; p += increment) { 00326 for (GiNaC::ex q=0; q <= one ; q += increment) { 00327 for (GiNaC::ex r=0; r <= one ; r += increment) { 00328 GiNaC::ex eq = polynom == GiNaC::numeric(0); 00329 equations.append(eq.subs(GiNaC::lst(x == p, y == q, z==r))); 00330 dofs.push_back(GiNaC::lst(p,q,r)); 00331 } 00332 } 00333 } 00334 00335 / * 00336 GiNaC::ex subs = lsolve(equations, variables); 00337 Nj = polynom.subs(subs); 00338 Ns.push_back(Nj); 00339 */ 00340 } 00341 00342 // invert the matrix: 00343 // GiNaC has a bit strange way to invert a matrix. 00344 // It solves the system AA^{-1} = Id. 00345 // It seems that this way is the only way to do 00346 // properly with the solve_algo::gauss flag. 00347 // 00348 00349 // std::cout <<"no variables "<<variables.nops()<<std::endl; 00350 // std::cout <<"no equations "<<equations.nops()<<std::endl; 00351 // print(equations); 00352 // print(variables); 00353 GiNaC::matrix b; GiNaC::matrix A; 00354 matrix_from_equations(equations, variables, A, b); 00355 00356 unsigned int ncols = A.cols(); 00357 GiNaC::matrix vars_sq(ncols, ncols); 00358 00359 // matrix of symbols 00360 for (unsigned r=0; r<ncols; ++r) 00361 for (unsigned c=0; c<ncols; ++c) 00362 vars_sq(r, c) = GiNaC::symbol(); 00363 00364 GiNaC::matrix id(ncols, ncols); 00365 00366 // identity 00367 const GiNaC::ex _ex1(1); 00368 for (unsigned i=0; i<ncols; ++i) 00369 id(i, i) = _ex1; 00370 00371 // invert the matrix 00372 GiNaC::matrix m_inv(ncols, ncols); 00373 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss); 00374 00375 for (unsigned int i=0; i<dofs.size(); i++) 00376 { 00377 b.let_op(i) = GiNaC::numeric(1); 00378 GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b)); 00379 00380 GiNaC::lst subs; 00381 for (unsigned int ii=0; ii<xx.nops(); ii++) 00382 { 00383 subs.append(variables.op(ii) == xx.op(ii)); 00384 } 00385 GiNaC::ex Nj = polynom.subs(subs); 00386 Ns.insert(Ns.end(), Nj); 00387 b.let_op(i) = GiNaC::numeric(0); 00388 } 00389 } 00390 00391 // ------------VectorLagrange --- 00392 00393 VectorLagrange::VectorLagrange() : StandardFE() 00394 { 00395 description = "VectorLagrange"; 00396 } 00397 00398 VectorLagrange::VectorLagrange(Polygon& p, unsigned int order, unsigned int size_) : StandardFE(p, order) 00399 { 00400 size = size_ < 0 ? nsd: size_; 00401 compute_basis_functions(); 00402 } 00403 00404 void VectorLagrange:: compute_basis_functions() 00405 { 00406 00407 // remove previously computed basis functions and dofs 00408 Ns.clear(); 00409 dofs.clear(); 00410 00411 if ( order < 1 ) 00412 { 00413 throw(std::logic_error("Lagrangian elements must be of order 1 or higher.")); 00414 } 00415 00416 if ( p == NULL ) 00417 { 00418 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed")); 00419 } 00420 00421 if ( size == 0) 00422 { 00423 throw(std::logic_error("You need to set the size of the vector before the basisfunctions can be computed")); 00424 } 00425 00426 Lagrange fe; 00427 fe.set_order(order); 00428 fe.set_polygon(*p); 00429 fe.compute_basis_functions(); 00430 GiNaC::lst zero_list; 00431 for (unsigned int s=1; s<= size ; s++) 00432 { 00433 zero_list.append(0); 00434 } 00435 00436 for (unsigned int s=0; s< size ; s++) 00437 { 00438 for (unsigned int i=0; i< fe.nbf() ; i++) 00439 { 00440 GiNaC::lst Nis = zero_list; 00441 Nis.let_op(s) = fe.N(i); 00442 GiNaC::ex Nmat = GiNaC::matrix(size,1,Nis); 00443 Ns.insert(Ns.end(), Nmat); 00444 00445 GiNaC::lst dof = GiNaC::lst(fe.dof(i), s) ; 00446 dofs.insert(dofs.end(), dof); 00447 } 00448 } 00449 00450 description = "Vector" + fe.str(); 00451 } 00452 00453 void VectorLagrange:: set_size(unsigned int size_) 00454 { 00455 size = size_; 00456 } 00457 00458 // ------------TensorLagrange --- 00459 00460 TensorLagrange::TensorLagrange() : StandardFE() 00461 { 00462 description = "TensorLagrange"; 00463 } 00464 00465 TensorLagrange::TensorLagrange(Polygon& p, unsigned int order, unsigned int size_) : StandardFE(p, order) 00466 { 00467 size = size_ < 0 ? nsd: size_; 00468 compute_basis_functions(); 00469 } 00470 00471 void TensorLagrange:: compute_basis_functions() 00472 { 00473 00474 // remove previously computed basis functions and dofs 00475 Ns.clear(); 00476 dofs.clear(); 00477 00478 if ( order < 1 ) 00479 { 00480 throw(std::logic_error("Lagrangian elements must be of order 1 or higher.")); 00481 } 00482 00483 if ( p == NULL ) 00484 { 00485 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed")); 00486 } 00487 00488 if ( size == 0) 00489 { 00490 throw(std::logic_error("You need to set the size of the vector before the basisfunctions can be computed")); 00491 } 00492 00493 Lagrange fe; 00494 fe.set_order(order); 00495 fe.set_polygon(*p); 00496 fe.compute_basis_functions(); 00497 GiNaC::lst zero_list; 00498 for (unsigned int s=1; s<= size*size ; s++) 00499 { 00500 zero_list.append(0); 00501 } 00502 00503 for (unsigned int r=0; r< size ; r++) 00504 { 00505 for (unsigned int s=0; s< size ; s++) 00506 { 00507 for (unsigned int i=0; i< fe.nbf() ; i++) 00508 { 00509 GiNaC::lst Nis = zero_list; 00510 Nis.let_op((size)*r + s) = fe.N(i); 00511 GiNaC::ex Nmat = GiNaC::matrix(size,size,Nis); 00512 Ns.insert(Ns.end(), Nmat); 00513 00514 GiNaC::lst dof = GiNaC::lst(fe.dof(i), r, s) ; 00515 dofs.insert(dofs.end(), dof); 00516 } 00517 } 00518 } 00519 00520 description = "Tensor" + fe.str(); 00521 } 00522 00523 void TensorLagrange:: set_size(unsigned int size_) 00524 { 00525 size = size_; 00526 } 00527 00528 GiNaC::ex lagrange(unsigned int order, Polygon& p, const std::string & a) 00529 { 00530 if ( order < 1 ) 00531 { 00532 throw(std::logic_error("Can not create polynomials of order less than 1!")); 00533 } 00534 00535 GiNaC::ex A; 00536 GiNaC::ex ret; 00537 GiNaC::lst basis; 00538 00539 Lagrange fe(p,order); 00540 A = get_symbolic_matrix(1, fe.nbf(), a); 00541 00542 for (unsigned int i=0; i<fe.nbf(); i++) 00543 { 00544 ret += A.op(i)*fe.N(i); 00545 basis.append(fe.N(i)); 00546 } 00547 return GiNaC::lst(ret,matrix_to_lst2(A),basis); 00548 } 00549 00550 GiNaC::lst lagrangev(unsigned int no_fields, unsigned int order, Polygon& p, const std::string & a) 00551 { 00552 if ( order < 1 ) 00553 { 00554 throw(std::logic_error("Can not create polynomials of order less than 1!")); 00555 } 00556 00557 GiNaC::ex A; 00558 GiNaC::ex ret; 00559 GiNaC::lst basis; 00560 00561 VectorLagrange fe(p,order); 00562 A = get_symbolic_matrix(1, fe.nbf(), a); 00563 00564 for (unsigned int i=0; i<fe.nbf(); i++) 00565 { 00566 ret += A.op(i)*fe.N(i); 00567 basis.append(fe.N(i)); 00568 } 00569 return GiNaC::lst(ret,matrix_to_lst2(A),basis); 00570 } 00571 00572 } // namespace SyFi
1.7.6.1 
