#include <Lagrange.h>
Public Member Functions | |
Lagrange () | |
Lagrange (Polygon &p, unsigned int order=1) | |
virtual | ~Lagrange () |
virtual void | compute_basis_functions () |
def | __init__ |
def | compute_basis_functions |
Public Attributes | |
this | |
Static Private Attributes | |
dictionary | __swig_setmethods__ = {} |
tuple | __setattr__ = lambdaself,name,value:_swig_setattr(self, Lagrange, name, value) |
dictionary | __swig_getmethods__ = {} |
tuple | __getattr__ = lambdaself,name:_swig_getattr(self, Lagrange, name) |
__repr__ = _swig_repr | |
__swig_destroy__ = _SyFi.delete_Lagrange | |
__del__ = lambdaself:None; |
Proxy of C++ SyFi::Lagrange class
Definition at line 12 of file Lagrange.h.
SyFi::Lagrange::Lagrange | ( | ) |
Definition at line 17 of file Lagrange.cpp.
References SyFi::StandardFE::description.
00017 : StandardFE() 00018 { 00019 description = "Lagrange"; 00020 }
SyFi::Lagrange::Lagrange | ( | Polygon & | p, | |
unsigned int | order = 1 | |||
) |
Definition at line 22 of file Lagrange.cpp.
References compute_basis_functions().
00022 : StandardFE(p, order) 00023 { 00024 compute_basis_functions(); 00025 }
virtual SyFi::Lagrange::~Lagrange | ( | ) | [inline, virtual] |
def SyFi::Lagrange::__init__ | ( | self, | ||
args | ||||
) |
__init__(self) -> Lagrange __init__(self, Polygon p, unsigned int order = 1) -> Lagrange __init__(self, Polygon p) -> Lagrange
Reimplemented from SyFi::StandardFE.
Reimplemented in SyFi::DiscontinuousLagrange.
Definition at line 2120 of file SyFi.py.
02120 : 02121 """ 02122 __init__(self) -> Lagrange 02123 __init__(self, Polygon p, unsigned int order = 1) -> Lagrange 02124 __init__(self, Polygon p) -> Lagrange 02125 """ 02126 this = _SyFi.new_Lagrange(*args) 02127 try: self.this.append(this) 02128 except: self.this = this
def SyFi::Lagrange::compute_basis_functions | ( | self | ) |
compute_basis_functions(self)
Reimplemented from SyFi::StandardFE.
Reimplemented in SyFi::DiscontinuousLagrange.
Definition at line 2131 of file SyFi.py.
02131 : 02132 """compute_basis_functions(self)""" 02133 return _SyFi.Lagrange_compute_basis_functions(self) 02134 Lagrange_swigregister = _SyFi.Lagrange_swigregister
void SyFi::Lagrange::compute_basis_functions | ( | ) | [virtual] |
Reimplemented from SyFi::StandardFE.
Reimplemented in SyFi::DiscontinuousLagrange.
Definition at line 27 of file Lagrange.cpp.
References SyFi::bernstein(), SyFi::bezier_ordinates(), SyFi::StandardFE::description, SyFi::StandardFE::dof(), SyFi::StandardFE::dofs, demos::crouzeixraviart::fe, SyFi::istr(), SyFi_polygons::spacetimedomain::l, SyFi::matrix_from_equations(), SyFi::StandardFE::N(), SyFi::StandardFE::nbf(), SyFi::StandardFE::Ns, SyFi::StandardFE::order, SyFi::StandardFE::p, SyFi::pol(), SyFi::Polygon::str(), SyFi::t, SyFi::x, SyFi::y, and SyFi::z.
Referenced by SyFi::TensorLagrange::compute_basis_functions(), SyFi::VectorLagrange::compute_basis_functions(), Lagrange(), and main().
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 }
SyFi::Lagrange::__del__ = lambdaself:None; [static, private] |
tuple SyFi::Lagrange::__getattr__ = lambdaself,name:_swig_getattr(self, Lagrange, name) [static, private] |
SyFi::Lagrange::__repr__ = _swig_repr [static, private] |
tuple SyFi::Lagrange::__setattr__ = lambdaself,name,value:_swig_setattr(self, Lagrange, name, value) [static, private] |
SyFi::Lagrange::__swig_destroy__ = _SyFi.delete_Lagrange [static, private] |
dictionary SyFi::Lagrange::__swig_getmethods__ = {} [static, private] |
dictionary SyFi::Lagrange::__swig_setmethods__ = {} [static, private] |