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