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