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 "Nedelec.h" 00005 #include <fstream> 00006 #include "tools.h" 00007 00008 using std::cout; 00009 using std::endl; 00010 using std::string; 00011 using GiNaC::exmap; 00012 00013 namespace SyFi 00014 { 00015 00016 Nedelec:: Nedelec() : StandardFE() 00017 { 00018 description = "Nedelec"; 00019 } 00020 00021 Nedelec:: Nedelec(Polygon& p, int order) : StandardFE(p, order) 00022 { 00023 compute_basis_functions(); 00024 } 00025 00026 void Nedelec:: compute_basis_functions() 00027 { 00028 00029 // remove previously computed basis functions and dofs 00030 Ns.clear(); 00031 dofs.clear(); 00032 00033 if ( order < 0 ) 00034 { 00035 throw(std::logic_error("Nedelec elements must be of order 0 or higher.")); 00036 } 00037 00038 if ( p == NULL ) 00039 { 00040 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed")); 00041 } 00042 00043 if ( p->str().find("Line") != string::npos ) 00044 { 00045 00046 cout <<"Can not define the Nedelec element on a line"<<endl; 00047 00048 } 00049 else if ( p->str().find("Triangle") != string::npos ) 00050 { 00051 00052 description = istr("Nedelec_", order) + "2D"; 00053 00054 int k = order; 00055 int removed_dofs=0; 00056 00057 Triangle& triangle = (Triangle&)(*p); 00058 GiNaC::lst equations; 00059 GiNaC::lst variables; 00060 00061 // create r 00062 GiNaC::ex R_k = homogenous_polv(2,k+1, 2, "a"); 00063 GiNaC::ex R_k_x = R_k.op(0).op(0); 00064 GiNaC::ex R_k_y = R_k.op(0).op(1); 00065 00066 // Equations that make sure that r*x = 0 00067 GiNaC::ex rx = (R_k_x*x + R_k_y*y).expand(); 00068 exmap pol_map = pol2basisandcoeff(rx); 00069 exmap::iterator iter; 00070 for (iter = pol_map.begin(); iter != pol_map.end(); iter++) 00071 { 00072 if ((*iter).second != 0 ) 00073 { 00074 equations.append((*iter).second == 0 ); 00075 removed_dofs++; 00076 } 00077 } 00078 00079 // create p 00080 GiNaC::ex P_k = bernsteinv(2,k, triangle, "b"); 00081 GiNaC::ex P_k_x = P_k.op(0).op(0); 00082 GiNaC::ex P_k_y = P_k.op(0).op(1); 00083 00084 // collect the variables of r and p in one list 00085 variables = collapse(GiNaC::lst(collapse(GiNaC::ex_to<GiNaC::lst>(R_k.op(1))), 00086 collapse(GiNaC::ex_to<GiNaC::lst>(P_k.op(1))))); 00087 00088 // create the polynomial space 00089 GiNaC::lst pspace = GiNaC::lst( R_k_x + P_k_x, 00090 R_k_y + P_k_y); 00091 00092 int counter = 0; 00093 GiNaC::symbol t("t"); 00094 GiNaC::ex dofi; 00095 // dofs related to edges 00096 for (int i=0; i< 3; i++) 00097 { 00098 Line line = triangle.line(i); 00099 GiNaC::lst tangent_vec = tangent(triangle, i); 00100 GiNaC::ex bernstein_pol = bernstein(order, line, istr("a",i)); 00101 GiNaC::ex basis_space = bernstein_pol.op(2); 00102 GiNaC::ex pspace_t = inner(pspace, tangent_vec); 00103 00104 GiNaC::ex basis; 00105 for (unsigned int j=0; j< basis_space.nops(); j++) 00106 { 00107 counter++; 00108 basis = basis_space.op(j); 00109 GiNaC::ex integrand = pspace_t*basis; 00110 dofi = line.integrate(integrand); 00111 GiNaC::ex eq = dofi == GiNaC::numeric(0); 00112 equations.append(eq); 00113 00114 GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),j); 00115 dofs.insert(dofs.end(), d); 00116 00117 } 00118 } 00119 00120 // dofs related to the whole triangle 00121 GiNaC::lst bernstein_polv; 00122 if ( order > 0) 00123 { 00124 counter++; 00125 bernstein_polv = bernsteinv(2,order-1, triangle, "a"); 00126 GiNaC::ex basis_space = bernstein_polv.op(2); 00127 for (unsigned int i=0; i< basis_space.nops(); i++) 00128 { 00129 GiNaC::lst basis = GiNaC::ex_to<GiNaC::lst>(basis_space.op(i)); 00130 GiNaC::ex integrand = inner(pspace, basis); 00131 dofi = triangle.integrate(integrand); 00132 GiNaC::ex eq = dofi == GiNaC::numeric(0); 00133 equations.append(eq); 00134 00135 GiNaC::lst d = GiNaC::lst(GiNaC::lst(triangle.vertex(0), triangle.vertex(1), triangle.vertex(2)), i); 00136 dofs.insert(dofs.end(), d); 00137 } 00138 } 00139 00140 // invert the matrix: 00141 // GiNaC has a bit strange way to invert a matrix. 00142 // It solves the system AA^{-1} = Id. 00143 // It seems that this way is the only way to do 00144 // properly with the solve_algo::gauss flag. 00145 GiNaC::matrix b; GiNaC::matrix A; 00146 matrix_from_equations(equations, variables, A, b); 00147 00148 unsigned int ncols = A.cols(); 00149 GiNaC::matrix vars_sq(ncols, ncols); 00150 00151 // matrix of symbols 00152 for (unsigned r=0; r<ncols; ++r) 00153 for (unsigned c=0; c<ncols; ++c) 00154 vars_sq(r, c) = GiNaC::symbol(); 00155 00156 GiNaC::matrix id(ncols, ncols); 00157 00158 // identity 00159 const GiNaC::ex _ex1(1); 00160 for (unsigned i=0; i<ncols; ++i) 00161 id(i, i) = _ex1; 00162 00163 // invert the matrix 00164 GiNaC::matrix m_inv(ncols, ncols); 00165 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss); 00166 00167 for (unsigned int i=0; i<dofs.size(); i++) 00168 { 00169 b.let_op(removed_dofs + i) = GiNaC::numeric(1); 00170 GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b)); 00171 00172 GiNaC::lst subs; 00173 for (unsigned int ii=0; ii<xx.nops(); ii++) 00174 { 00175 subs.append(variables.op(ii) == xx.op(ii)); 00176 } 00177 GiNaC::ex Nj1 = pspace.op(0).subs(subs); 00178 GiNaC::ex Nj2 = pspace.op(1).subs(subs); 00179 Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2))); 00180 b.let_op(removed_dofs + i) = GiNaC::numeric(0); 00181 } 00182 00183 } 00184 else if ( p->str().find("Tetrahedron") != string::npos ) 00185 { 00186 00187 description = istr("Nedelec_", order) + "3D"; 00188 00189 int k = order; 00190 int removed_dofs=0; 00191 00192 Tetrahedron& tetrahedron= (Tetrahedron&)(*p); 00193 GiNaC::lst equations; 00194 GiNaC::lst variables; 00195 00196 // create r 00197 GiNaC::ex R_k = homogenous_polv(3,k+1, 3, "a"); 00198 GiNaC::ex R_k_x = R_k.op(0).op(0); 00199 GiNaC::ex R_k_y = R_k.op(0).op(1); 00200 GiNaC::ex R_k_z = R_k.op(0).op(2); 00201 00202 // Equations that make sure that r*x = 0 00203 GiNaC::ex rx = (R_k_x*x + R_k_y*y + R_k_z*z).expand(); 00204 exmap pol_map = pol2basisandcoeff(rx); 00205 exmap::iterator iter; 00206 for (iter = pol_map.begin(); iter != pol_map.end(); iter++) 00207 { 00208 if ((*iter).second != 0 ) 00209 { 00210 equations.append((*iter).second == 0 ); 00211 removed_dofs++; 00212 } 00213 } 00214 00215 // create p 00216 GiNaC::ex P_k = bernsteinv(3,k, tetrahedron, "b"); 00217 GiNaC::ex P_k_x = P_k.op(0).op(0); 00218 GiNaC::ex P_k_y = P_k.op(0).op(1); 00219 GiNaC::ex P_k_z = P_k.op(0).op(2); 00220 00221 // collect the variables of r and p in one list 00222 variables = collapse(GiNaC::lst(collapse(GiNaC::ex_to<GiNaC::lst>(R_k.op(1))), 00223 collapse(GiNaC::ex_to<GiNaC::lst>(P_k.op(1))))); 00224 00225 // create the polynomial space 00226 GiNaC::lst pspace = GiNaC::lst( R_k_x + P_k_x, 00227 R_k_y + P_k_y, 00228 R_k_z + P_k_z); 00229 00230 int counter = 0; 00231 GiNaC::symbol t("t"); 00232 GiNaC::ex dofi; 00233 00234 // dofs related to edges 00235 for (int i=0; i< 6; i++) 00236 { 00237 Line line = tetrahedron.line(i); 00238 GiNaC::ex line_repr = line.repr(t); 00239 GiNaC::lst tangent_vec = GiNaC::lst(line_repr.op(0).rhs().coeff(t,1), 00240 line_repr.op(1).rhs().coeff(t,1), 00241 line_repr.op(2).rhs().coeff(t,1)); 00242 00243 GiNaC::ex bernstein_pol = bernstein(order, line, istr("a",i)); 00244 GiNaC::ex basis_space = bernstein_pol.op(2); 00245 GiNaC::ex pspace_t = inner(pspace, tangent_vec); 00246 00247 GiNaC::ex basis; 00248 for (unsigned int j=0; j< basis_space.nops(); j++) 00249 { 00250 counter++; 00251 basis = basis_space.op(j); 00252 GiNaC::ex integrand = pspace_t*basis; 00253 dofi = line.integrate(integrand); 00254 GiNaC::ex eq = dofi == GiNaC::numeric(0); 00255 equations.append(eq); 00256 00257 GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),j); 00258 dofs.insert(dofs.end(), d); 00259 00260 } 00261 } 00262 00263 // dofs related to faces 00264 if ( order > 0 ) 00265 { 00266 for (int i=0; i< 4; i++) 00267 { 00268 Triangle triangle = tetrahedron.triangle(i); 00269 GiNaC::ex bernstein_pol = bernsteinv(3,order-1, triangle, istr("b", i)); 00270 GiNaC::ex basis_space = bernstein_pol.op(2); 00271 00272 GiNaC::ex basis; 00273 GiNaC::lst normal_vec = normal(tetrahedron, i); 00274 GiNaC::ex pspace_n = cross(pspace, normal_vec); 00275 if ( normal_vec.op(0) != GiNaC::numeric(0) && 00276 normal_vec.op(1) != GiNaC::numeric(0) && 00277 normal_vec.op(2) != GiNaC::numeric(0) ) 00278 { 00279 for (unsigned int j=0; j<basis_space.nops(); j++) 00280 { 00281 basis = basis_space.op(j); 00282 if ( basis.op(0) != 0 || basis.op(1) != 0 ) 00283 { 00284 GiNaC::ex integrand = inner(pspace_n,basis); 00285 if ( integrand != GiNaC::numeric(0) ) 00286 { 00287 dofi = triangle.integrate(integrand); 00288 GiNaC::ex eq = dofi == GiNaC::numeric(0); 00289 equations.append(eq); 00290 } 00291 } 00292 } 00293 00294 } 00295 else 00296 { 00297 for (unsigned int j=0; j<basis_space.nops(); j++) 00298 { 00299 basis = basis_space.op(j); 00300 GiNaC::ex integrand = inner(pspace_n,basis); 00301 if ( integrand != GiNaC::numeric(0) ) 00302 { 00303 dofi = triangle.integrate(integrand); 00304 GiNaC::ex eq = dofi == GiNaC::numeric(0); 00305 equations.append(eq); 00306 } 00307 } 00308 } 00309 } 00310 } 00311 00312 // dofs related to tetrahedron 00313 if ( order > 1 ) 00314 { 00315 GiNaC::ex bernstein_pol = bernsteinv(3,order-2, tetrahedron, istr("c", 0)); 00316 GiNaC::ex basis_space = bernstein_pol.op(2); 00317 GiNaC::ex basis; 00318 for (unsigned int j=0; j<basis_space.nops(); j++) 00319 { 00320 basis = basis_space.op(j); 00321 GiNaC::ex integrand = inner(pspace,basis); 00322 dofi = tetrahedron.integrate(integrand); 00323 GiNaC::ex eq = dofi == GiNaC::numeric(0); 00324 equations.append(eq); 00325 00326 GiNaC::lst d = GiNaC::lst(GiNaC::lst(tetrahedron.vertex(0), tetrahedron.vertex(1), tetrahedron.vertex(2), tetrahedron.vertex(3)), j); 00327 dofs.insert(dofs.end(), d); 00328 00329 } 00330 } 00331 00332 // invert the matrix: 00333 // GiNaC has a bit strange way to invert a matrix. 00334 // It solves the system AA^{-1} = Id. 00335 // It seems that this way is the only way to do 00336 // properly with the solve_algo::gauss flag. 00337 GiNaC::matrix b; GiNaC::matrix A; 00338 matrix_from_equations(equations, variables, A, b); 00339 00340 unsigned int ncols = A.cols(); 00341 GiNaC::matrix vars_sq(ncols, ncols); 00342 00343 // matrix of symbols 00344 for (unsigned r=0; r<ncols; ++r) 00345 for (unsigned c=0; c<ncols; ++c) 00346 vars_sq(r, c) = GiNaC::symbol(); 00347 00348 GiNaC::matrix id(ncols, ncols); 00349 00350 // identity 00351 const GiNaC::ex _ex1(1); 00352 for (unsigned i=0; i<ncols; ++i) 00353 id(i, i) = _ex1; 00354 00355 // invert the matrix 00356 GiNaC::matrix m_inv(ncols, ncols); 00357 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss); 00358 00359 for (unsigned int i=0; i<dofs.size(); i++) 00360 { 00361 b.let_op(removed_dofs + i) = GiNaC::numeric(1); 00362 GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b)); 00363 00364 GiNaC::lst subs; 00365 for (unsigned int ii=0; ii<xx.nops(); ii++) 00366 { 00367 subs.append(variables.op(ii) == xx.op(ii)); 00368 } 00369 GiNaC::ex Nj1 = pspace.op(0).subs(subs); 00370 GiNaC::ex Nj2 = pspace.op(1).subs(subs); 00371 GiNaC::ex Nj3 = pspace.op(2).subs(subs); 00372 Ns.insert(Ns.end(), GiNaC::matrix(3,1,GiNaC::lst(Nj1,Nj2,Nj3))); 00373 b.let_op(removed_dofs + i) = GiNaC::numeric(0); 00374 } 00375 00376 } 00377 00378 } 00379 00380 } // namespace SyFi