#include <Nedelec2Hdiv.h>
Public Member Functions | |
Nedelec2Hdiv () | |
Nedelec2Hdiv (Polygon &p, unsigned int order=1) | |
virtual | ~Nedelec2Hdiv () |
virtual void | compute_basis_functions () |
def | __init__ |
def | compute_basis_functions |
Public Attributes | |
GiNaC::lst | dof_repr |
this | |
Static Private Attributes | |
dictionary | __swig_setmethods__ = {} |
tuple | __setattr__ = lambdaself,name,value:_swig_setattr(self, Nedelec2Hdiv, name, value) |
dictionary | __swig_getmethods__ = {} |
tuple | __getattr__ = lambdaself,name:_swig_getattr(self, Nedelec2Hdiv, name) |
__repr__ = _swig_repr | |
__swig_destroy__ = _SyFi.delete_Nedelec2Hdiv | |
__del__ = lambdaself:None; |
Proxy of C++ SyFi::Nedelec2Hdiv class
Definition at line 12 of file Nedelec2Hdiv.h.
SyFi::Nedelec2Hdiv::Nedelec2Hdiv | ( | ) |
Definition at line 18 of file Nedelec2Hdiv.cpp.
References SyFi::StandardFE::description.
00018 : StandardFE() 00019 { 00020 description = "Nedelec2Hdiv"; 00021 }
SyFi::Nedelec2Hdiv::Nedelec2Hdiv | ( | Polygon & | p, | |
unsigned int | order = 1 | |||
) |
Definition at line 23 of file Nedelec2Hdiv.cpp.
References compute_basis_functions().
00023 : StandardFE(p, order) 00024 { 00025 compute_basis_functions(); 00026 }
virtual SyFi::Nedelec2Hdiv::~Nedelec2Hdiv | ( | ) | [inline, virtual] |
def SyFi::Nedelec2Hdiv::__init__ | ( | self, | ||
args | ||||
) |
__init__(self) -> Nedelec2Hdiv __init__(self, Polygon p, unsigned int order = 1) -> Nedelec2Hdiv __init__(self, Polygon p) -> Nedelec2Hdiv
Reimplemented from SyFi::StandardFE.
Definition at line 2582 of file SyFi.py.
02582 : 02583 """ 02584 __init__(self) -> Nedelec2Hdiv 02585 __init__(self, Polygon p, unsigned int order = 1) -> Nedelec2Hdiv 02586 __init__(self, Polygon p) -> Nedelec2Hdiv 02587 """ 02588 this = _SyFi.new_Nedelec2Hdiv(*args) 02589 try: self.this.append(this) 02590 except: self.this = this
def SyFi::Nedelec2Hdiv::compute_basis_functions | ( | self | ) |
compute_basis_functions(self)
Reimplemented from SyFi::StandardFE.
Definition at line 2593 of file SyFi.py.
02593 : 02594 """compute_basis_functions(self)""" 02595 return _SyFi.Nedelec2Hdiv_compute_basis_functions(self) 02596 Nedelec2Hdiv_swigregister = _SyFi.Nedelec2Hdiv_swigregister
void SyFi::Nedelec2Hdiv::compute_basis_functions | ( | ) | [virtual] |
Reimplemented from SyFi::StandardFE.
Definition at line 28 of file Nedelec2Hdiv.cpp.
References SyFi::exmap::begin(), SyFi::bernstein(), SyFi::bernsteinv(), SyFi::coeff(), SyFi::collapse(), SyFi::StandardFE::description, dof_repr, SyFi::StandardFE::dofs, SyFi::exmap::end(), SyFi::homogenous_polv(), SyFi::inner(), demos::poisson1::integrand, SyFi::Tetrahedron::integrate(), SyFi::Triangle::integrate(), SyFi::istr(), SyFi::exmap::iterator(), SyFi::matrix_from_equations(), test::n, SyFi::normal(), SyFi::StandardFE::Ns, SyFi::StandardFE::order, SyFi::StandardFE::p, SyFi::pol2basisandcoeff(), SyFi::Polygon::str(), SyFi::t, demos::geom_test::tetrahedron, SyFi::Tetrahedron::triangle(), demos::crouzeixraviart::triangle, fem_sympy::u, SyFi::Polygon::vertex(), SyFi::x, SyFi::y, and SyFi::z.
Referenced by SyFi::ArnoldFalkWintherWeakSymSigma::compute_basis_functions(), main(), and Nedelec2Hdiv().
00029 { 00030 00031 // remove previously computed basis functions and dofs 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 // create p 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 // dofs related to edges 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 // GiNaC::lst d = GiNaC::lst(triangle.integrate(x*basis), 00102 // triangle.integrate(y*basis), 00103 // triangle.integrate(z*basis)); 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 // dofs related to tetrahedron 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 // GiNaC::lst d = GiNaC::lst(tetrahedron.integrate(x*basis.op(0)), 00137 // tetrahedron.integrate(y*basis.op(1)), 00138 // tetrahedron.integrate(z*basis.op(2))); 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 // Construction of S_k 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 // Equations that make sure that r*x = 0 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 // 2 basis functions is added 00193 // The first: 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 // The second: 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 // dofs related to tetrahedron 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 // invert the matrix: 00239 // GiNaC has a bit strange way to invert a matrix. 00240 // It solves the system AA^{-1} = Id. 00241 // It seems that this way is the only way to do 00242 // properly with the solve_algo::gauss flag. 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 // matrix of symbols 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 // identity 00258 const GiNaC::ex _ex1(1); 00259 for (unsigned i=0; i<ncols; ++i) 00260 id(i, i) = _ex1; 00261 00262 // invert the matrix 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 }
SyFi::Nedelec2Hdiv::__del__ = lambdaself:None; [static, private] |
tuple SyFi::Nedelec2Hdiv::__getattr__ = lambdaself,name:_swig_getattr(self, Nedelec2Hdiv, name) [static, private] |
SyFi::Nedelec2Hdiv::__repr__ = _swig_repr [static, private] |
tuple SyFi::Nedelec2Hdiv::__setattr__ = lambdaself,name,value:_swig_setattr(self, Nedelec2Hdiv, name, value) [static, private] |
SyFi::Nedelec2Hdiv::__swig_destroy__ = _SyFi.delete_Nedelec2Hdiv [static, private] |
dictionary SyFi::Nedelec2Hdiv::__swig_getmethods__ = {} [static, private] |
dictionary SyFi::Nedelec2Hdiv::__swig_setmethods__ = {} [static, private] |
GiNaC::lst SyFi::Nedelec2Hdiv::dof_repr |