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