00001
00002
00003
00004 #include "Nedelec2Hdiv.h"
00005 #include <fstream>
00006 #include "tools.h"
00007
00008 using std::cout;
00009 using std::endl;
00010 using std::string;
00011
00012 using GiNaC::exmap;
00013 using namespace GiNaC;
00014
00015 namespace SyFi
00016 {
00017
00018 Nedelec2Hdiv:: Nedelec2Hdiv() : StandardFE()
00019 {
00020 description = "Nedelec2Hdiv";
00021 }
00022
00023 Nedelec2Hdiv:: Nedelec2Hdiv(Polygon& p, unsigned int order) : StandardFE(p, order)
00024 {
00025 compute_basis_functions();
00026 }
00027
00028 void Nedelec2Hdiv:: compute_basis_functions()
00029 {
00030
00031
00032 Ns.clear();
00033 dofs.clear();
00034
00035 if ( order < 1 )
00036 {
00037 throw(std::logic_error("Nedelec2Hdiv elements must be of order 1 or higher."));
00038 }
00039
00040 if ( p == NULL )
00041 {
00042 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00043 }
00044
00045 if ( p->str().find("Line") != string::npos )
00046 {
00047
00048 cout <<"Can not define the Nedelec2Hdiv element on a line"<<endl;
00049
00050 }
00051 else if ( p->str().find("Triangle") != string::npos )
00052 {
00053
00054 cout <<"Can not define the Nedelec2Hdiv element on a Triangle "<<endl;
00055
00056 }
00057 else if ( p->str().find("Tetrahedron") != string::npos )
00058 {
00059
00060 description = istr( "Nedelec2Hdiv_", order) + "_3D";
00061
00062 int k = order;
00063
00064 Tetrahedron& tetrahedron= (Tetrahedron&)(*p);
00065 GiNaC::lst equations;
00066 GiNaC::lst variables;
00067
00068
00069 GiNaC::ex P_k = bernsteinv(3,k, tetrahedron, "b");
00070 GiNaC::ex P_k_x = P_k.op(0).op(0);
00071 GiNaC::ex P_k_y = P_k.op(0).op(1);
00072 GiNaC::ex P_k_z = P_k.op(0).op(2);
00073
00074 GiNaC::lst pspace = GiNaC::lst( P_k_x, P_k_y, P_k_z);
00075
00076 variables = collapse(GiNaC::ex_to<GiNaC::lst>(P_k.op(1)));
00077
00078 int counter = 0;
00079 GiNaC::symbol t("t");
00080 GiNaC::ex dofi;
00081 GiNaC::ex bernstein_pol;
00082
00083
00084 for (int i=0; i< 4; i++)
00085 {
00086 Triangle triangle = tetrahedron.triangle(i);
00087 GiNaC::lst normal_vec = normal(tetrahedron, i);
00088 bernstein_pol = bernstein(order, triangle, istr("a",i));
00089 GiNaC::ex basis_space = bernstein_pol.op(2);
00090 GiNaC::ex pspace_n = inner(pspace, normal_vec);
00091
00092 GiNaC::ex basis;
00093 for (unsigned int j=0; j< basis_space.nops(); j++)
00094 {
00095 counter++;
00096 basis = basis_space.op(j);
00097 GiNaC::ex integrand = pspace_n*basis;
00098 dofi = triangle.integrate(integrand);
00099 GiNaC::ex eq = dofi == GiNaC::numeric(0);
00100 equations.append(eq);
00101
00102
00103
00104 GiNaC::lst d = GiNaC::lst(GiNaC::lst(triangle.vertex(0),
00105 triangle.vertex(1),
00106 triangle.vertex(2)),j);
00107
00108 dofs.insert(dofs.end(), d);
00109
00110 GiNaC::ex u = GiNaC::matrix(3,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]"), GiNaC::symbol("v[2]")));
00111 GiNaC::ex n = GiNaC::matrix(3,1,GiNaC::lst(GiNaC::symbol("normal_vec[0]"), GiNaC::symbol("normal_vec[1]"), GiNaC::symbol("normal_vec[2]")));
00112 dof_repr.append(GiNaC::lst(inner(u,n)*basis.subs( x == GiNaC::symbol("xi[0]"))
00113 .subs( y == GiNaC::symbol("xi[1]"))
00114 .subs( z == GiNaC::symbol("xi[2]")), d));
00115
00116 }
00117 }
00118
00119
00120
00121 int tetradofs = 0;
00122 if ( order > 1 )
00123 {
00124 GiNaC::ex bernstein_pol = bernsteinv(3,k-2, tetrahedron, istr("c", 0));
00125 GiNaC::ex basis_space = bernstein_pol.op(2);
00126 GiNaC::ex basis;
00127 tetradofs += basis_space.nops();
00128 for (unsigned int j=0; j<basis_space.nops(); j++)
00129 {
00130 basis = basis_space.op(j);
00131 GiNaC::ex integrand = inner(pspace,basis);
00132 dofi = tetrahedron.integrate(integrand);
00133 GiNaC::ex eq = dofi == GiNaC::numeric(0);
00134 equations.append(eq);
00135
00136
00137
00138
00139
00140 GiNaC::lst d = GiNaC::lst(tetrahedron.vertex(0),
00141 tetrahedron.vertex(1),
00142 tetrahedron.vertex(2),
00143 tetrahedron.vertex(3),j);
00144
00145 dofs.insert(dofs.end(), d);
00146
00147 }
00148 }
00149
00150
00151
00152
00153 if ( order >= 1 )
00154 {
00155 GiNaC::ex H_k = homogenous_polv(3,k-1, 3, "a");
00156 GiNaC::ex H_k_x = H_k.op(0).op(0);
00157 GiNaC::ex H_k_y = H_k.op(0).op(1);
00158 GiNaC::ex H_k_z = H_k.op(0).op(2);
00159
00160 GiNaC::lst H_variables = collapse(GiNaC::ex_to<GiNaC::lst>(H_k.op(1)));
00161
00162
00163
00164 GiNaC::ex rx = (H_k_x*x + H_k_y*y + H_k_z*z).expand();
00165 exmap pol_map = pol2basisandcoeff(rx);
00166 exmap::iterator iter;
00167 GiNaC::lst S_k;
00168 GiNaC::lst S_k_equations;
00169
00170 GiNaC::lst null_eqs;
00171 for (unsigned int i=0; i<H_variables.nops(); i++)
00172 {
00173 null_eqs.append( H_variables.op(i) == 0);
00174 }
00175
00176 for (iter = pol_map.begin(); iter != pol_map.end(); iter++)
00177 {
00178 GiNaC::ex coeff = (*iter).second;
00179 GiNaC::ex basis;
00180 if (coeff.nops() > 1 )
00181 {
00182 if (coeff.nops() == 2)
00183 {
00184 S_k_equations.remove_all();
00185 S_k_equations.append(coeff.op(0) == GiNaC::numeric(1));
00186 S_k_equations.append(coeff.op(1) == GiNaC::numeric(-1));
00187 basis = H_k.op(0).subs(S_k_equations).subs(null_eqs);;
00188 S_k.append(basis);
00189 }
00190 else if ( coeff.nops() == 3 )
00191 {
00192
00193
00194
00195 S_k_equations.remove_all();
00196 S_k_equations.append(coeff.op(0) == GiNaC::numeric(-1,2));
00197 S_k_equations.append(coeff.op(1) == GiNaC::numeric(1));
00198 S_k_equations.append(coeff.op(2) == GiNaC::numeric(-1,2));
00199 basis = H_k.op(0).subs(S_k_equations).subs(null_eqs);;
00200 S_k.append(basis);
00201
00202
00203 S_k_equations.remove_all();
00204 S_k_equations.append(coeff.op(0) == GiNaC::numeric(-1,2));
00205 S_k_equations.append(coeff.op(1) == GiNaC::numeric(-1,2));
00206 S_k_equations.append(coeff.op(2) == GiNaC::numeric(1));
00207 basis = H_k.op(0).subs(S_k_equations).subs(null_eqs);;
00208 S_k.append(basis);
00209 }
00210 }
00211 }
00212
00213 std::cout <<"len (S_k) " <<S_k.nops()<<std::endl;
00214
00215
00216 if ( order >= 1 )
00217 {
00218 GiNaC::ex basis;
00219 for (unsigned int j=0; j<S_k.nops(); j++)
00220 {
00221 basis = S_k.op(j);
00222 GiNaC::ex integrand = inner(pspace,basis);
00223 dofi = tetrahedron.integrate(integrand);
00224 GiNaC::ex eq = dofi == GiNaC::numeric(0);
00225 equations.append(eq);
00226
00227 GiNaC::lst d = GiNaC::lst(tetrahedron.vertex(0),
00228 tetrahedron.vertex(1),
00229 tetrahedron.vertex(2),
00230 tetrahedron.vertex(3), tetradofs + j);
00231
00232 dofs.insert(dofs.end(), d);
00233
00234 }
00235 }
00236 }
00237
00238
00239
00240
00241
00242
00243
00244 GiNaC::matrix b; GiNaC::matrix A;
00245 matrix_from_equations(equations, variables, A, b);
00246
00247 unsigned int ncols = A.cols();
00248 GiNaC::matrix vars_sq(ncols, ncols);
00249
00250
00251 for (unsigned r=0; r<ncols; ++r)
00252 for (unsigned c=0; c<ncols; ++c)
00253 vars_sq(r, c) = GiNaC::symbol();
00254
00255 GiNaC::matrix id(ncols, ncols);
00256
00257
00258 const GiNaC::ex _ex1(1);
00259 for (unsigned i=0; i<ncols; ++i)
00260 id(i, i) = _ex1;
00261
00262
00263 GiNaC::matrix m_inv(ncols, ncols);
00264 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss);
00265
00266 for (unsigned int i=0; i<dofs.size(); i++)
00267 {
00268 b.let_op(i) = GiNaC::numeric(1);
00269 GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b));
00270
00271 GiNaC::lst subs;
00272 for (unsigned int ii=0; ii<xx.nops(); ii++)
00273 {
00274 subs.append(variables.op(ii) == xx.op(ii));
00275 }
00276 GiNaC::ex Nj1 = pspace.op(0).subs(subs);
00277 GiNaC::ex Nj2 = pspace.op(1).subs(subs);
00278 GiNaC::ex Nj3 = pspace.op(2).subs(subs);
00279 Ns.insert(Ns.end(), GiNaC::matrix(3,1,GiNaC::lst(Nj1,Nj2,Nj3)));
00280 b.let_op(i) = GiNaC::numeric(0);
00281 }
00282 }
00283 }
00284
00285 }