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