00001
00002
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
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
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
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
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
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
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
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
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
00141
00142
00143
00144
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
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
00159 const GiNaC::ex _ex1(1);
00160 for (unsigned i=0; i<ncols; ++i)
00161 id(i, i) = _ex1;
00162
00163
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
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
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
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
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
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
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
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
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
00333
00334
00335
00336
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
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
00351 const GiNaC::ex _ex1(1);
00352 for (unsigned i=0; i<ncols; ++i)
00353 id(i, i) = _ex1;
00354
00355
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 }