SyFi::Nedelec Class Reference

#include <Nedelec.h>

Inheritance diagram for SyFi::Nedelec:

SyFi::StandardFE SyFi::StandardFE SyFi::FE SyFi::FE SyFi::FE SyFi::FE SyFi::_object SyFi::_object SyFi::_object SyFi::_object

List of all members.

Public Member Functions

 Nedelec ()
 Nedelec (Polygon &p, int order=1)
virtual ~Nedelec ()
virtual void compute_basis_functions ()
def __init__
def compute_basis_functions

Public Attributes

 this

Static Private Attributes

dictionary __swig_setmethods__ = {}
tuple __setattr__ = lambdaself,name,value:_swig_setattr(self, Nedelec, name, value)
dictionary __swig_getmethods__ = {}
tuple __getattr__ = lambdaself,name:_swig_getattr(self, Nedelec, name)
 __repr__ = _swig_repr
 __swig_destroy__ = _SyFi.delete_Nedelec
 __del__ = lambdaself:None;


Detailed Description

Proxy of C++ SyFi::Nedelec class

Definition at line 12 of file Nedelec.h.


Constructor & Destructor Documentation

SyFi::Nedelec::Nedelec (  ) 

Definition at line 16 of file Nedelec.cpp.

References SyFi::StandardFE::description.

00016                           : StandardFE()
00017         {
00018                 description = "Nedelec";
00019         }

SyFi::Nedelec::Nedelec ( Polygon p,
int  order = 1 
)

Definition at line 21 of file Nedelec.cpp.

References compute_basis_functions().

00021                                                : StandardFE(p, order)
00022         {
00023                 compute_basis_functions();
00024         }

virtual SyFi::Nedelec::~Nedelec (  )  [inline, virtual]

Definition at line 17 of file Nedelec.h.

00017 {}


Member Function Documentation

def SyFi::Nedelec::__init__ (   self,
  args 
)

__init__(self) -> Nedelec
__init__(self, Polygon p, int order = 1) -> Nedelec
__init__(self, Polygon p) -> Nedelec

Reimplemented from SyFi::StandardFE.

Definition at line 2552 of file SyFi.py.

02552                              : 
02553         """
02554         __init__(self) -> Nedelec
02555         __init__(self, Polygon p, int order = 1) -> Nedelec
02556         __init__(self, Polygon p) -> Nedelec
02557         """
02558         this = _SyFi.new_Nedelec(*args)
02559         try: self.this.append(this)
02560         except: self.this = this

def SyFi::Nedelec::compute_basis_functions (   self  ) 

compute_basis_functions(self)

Reimplemented from SyFi::StandardFE.

Definition at line 2563 of file SyFi.py.

02563                                      :
02564         """compute_basis_functions(self)"""
02565         return _SyFi.Nedelec_compute_basis_functions(self)
02566 
Nedelec_swigregister = _SyFi.Nedelec_swigregister

void SyFi::Nedelec::compute_basis_functions (  )  [virtual]

Reimplemented from SyFi::StandardFE.

Definition at line 26 of file Nedelec.cpp.

References SyFi::exmap::begin(), SyFi::bernstein(), SyFi::bernsteinv(), SyFi::collapse(), SyFi::cross(), SyFi::StandardFE::description, SyFi::StandardFE::dofs, SyFi::exmap::end(), SyFi::homogenous_polv(), SyFi::inner(), demos::poisson1::integrand, SyFi::Tetrahedron::integrate(), SyFi::Triangle::integrate(), SyFi::Line::integrate(), SyFi::istr(), SyFi::exmap::iterator(), SyFi::Tetrahedron::line(), SyFi::Triangle::line(), SyFi::matrix_from_equations(), SyFi::normal(), SyFi::StandardFE::Ns, SyFi::StandardFE::order, SyFi::StandardFE::p, SyFi::pol2basisandcoeff(), SyFi::Line::repr(), SyFi::Polygon::str(), SyFi::t, SyFi::tangent(), demos::geom_test::tetrahedron, SyFi::Tetrahedron::triangle(), demos::crouzeixraviart::triangle, SyFi::Polygon::vertex(), SyFi::x, SyFi::y, and SyFi::z.

Referenced by main(), and Nedelec().

00027         {
00028 
00029                 // remove previously computed basis functions and dofs
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                         // create r
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                         // Equations that make sure that r*x = 0
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                         // create p
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                         // collect the variables of r and p  in one list
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                         // create the polynomial space
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                         // dofs related to edges
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                         // dofs related to the whole triangle
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                         // invert the matrix:
00141                         // GiNaC has a bit strange way to invert a matrix.
00142                         // It solves the system AA^{-1} = Id.
00143                         // It seems that this way is the only way to do
00144                         // properly with the solve_algo::gauss flag.
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                         // matrix of symbols
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                         // identity
00159                         const GiNaC::ex _ex1(1);
00160                         for (unsigned i=0; i<ncols; ++i)
00161                                 id(i, i) = _ex1;
00162 
00163                         // invert the matrix
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                         // create r
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                         // Equations that make sure that r*x = 0
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                         // create p
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                         // collect the variables of r and p  in one list
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                         // create the polynomial space
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                         // dofs related to edges
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                         // dofs related to faces
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                         // dofs related to tetrahedron
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                         // invert the matrix:
00333                         // GiNaC has a bit strange way to invert a matrix.
00334                         // It solves the system AA^{-1} = Id.
00335                         // It seems that this way is the only way to do
00336                         // properly with the solve_algo::gauss flag.
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                         // matrix of symbols
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                         // identity
00351                         const GiNaC::ex _ex1(1);
00352                         for (unsigned i=0; i<ncols; ++i)
00353                                 id(i, i) = _ex1;
00354 
00355                         // invert the matrix
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         }


Member Data Documentation

SyFi::Nedelec::__del__ = lambdaself:None; [static, private]

Reimplemented from SyFi::StandardFE.

Definition at line 2562 of file SyFi.py.

tuple SyFi::Nedelec::__getattr__ = lambdaself,name:_swig_getattr(self, Nedelec, name) [static, private]

Reimplemented from SyFi::StandardFE.

Definition at line 2550 of file SyFi.py.

SyFi::Nedelec::__repr__ = _swig_repr [static, private]

Reimplemented from SyFi::StandardFE.

Definition at line 2551 of file SyFi.py.

tuple SyFi::Nedelec::__setattr__ = lambdaself,name,value:_swig_setattr(self, Nedelec, name, value) [static, private]

Reimplemented from SyFi::StandardFE.

Definition at line 2547 of file SyFi.py.

SyFi::Nedelec::__swig_destroy__ = _SyFi.delete_Nedelec [static, private]

Reimplemented from SyFi::StandardFE.

Definition at line 2561 of file SyFi.py.

dictionary SyFi::Nedelec::__swig_getmethods__ = {} [static, private]

Reimplemented from SyFi::StandardFE.

Definition at line 2548 of file SyFi.py.

dictionary SyFi::Nedelec::__swig_setmethods__ = {} [static, private]

Reimplemented from SyFi::StandardFE.

Definition at line 2545 of file SyFi.py.

Reimplemented from SyFi::StandardFE.

Definition at line 2560 of file SyFi.py.


The documentation for this class was generated from the following files:

Generated on Mon Aug 31 16:17:07 2009 for SyFi by  doxygen 1.5.9