00001
00002
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
00031
00032
00033
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
00057
00058
00059
00060 GiNaC::ex polynom_space = bernstein(order, *p, "a");
00061
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
00080
00081
00082
00083 GiNaC::ex polynom_space = bernstein(order, *p, "a");
00084
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
00112
00113
00114
00115 GiNaC::ex polynom_space = bernstein(order, *p, "b");
00116
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
00140
00141
00142 GiNaC::ex polynom_space = bernstein(order, *p, "b");
00143
00144
00145
00146
00147
00148
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
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
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
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213 }
00214 else if ( p->str().find("ReferenceTetrahedron") != string::npos)
00215 {
00216
00217 description = istr("Lagrange_", order) + "_3D";
00218
00219
00220
00221
00222
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
00255
00256 GiNaC::ex polynom_space = pol(order, 3, "a");
00257
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
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326 }
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
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
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
00353 const GiNaC::ex _ex1(1);
00354 for (unsigned i=0; i<ncols; ++i)
00355 id(i, i) = _ex1;
00356
00357
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
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
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
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
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 }