SyFi::Nedelec2Hdiv Class Reference

#include <Nedelec2Hdiv.h>

Inheritance diagram for SyFi::Nedelec2Hdiv:

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

 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;


Detailed Description

Proxy of C++ SyFi::Nedelec2Hdiv class

Definition at line 12 of file Nedelec2Hdiv.h.


Constructor & Destructor Documentation

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]

Definition at line 18 of file Nedelec2Hdiv.h.

00018 {}


Member Function Documentation

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         }


Member Data Documentation

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

Reimplemented from SyFi::StandardFE.

Definition at line 2592 of file SyFi.py.

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

Reimplemented from SyFi::StandardFE.

Definition at line 2577 of file SyFi.py.

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

Reimplemented from SyFi::StandardFE.

Definition at line 2578 of file SyFi.py.

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

Reimplemented from SyFi::StandardFE.

Definition at line 2574 of file SyFi.py.

SyFi::Nedelec2Hdiv::__swig_destroy__ = _SyFi.delete_Nedelec2Hdiv [static, private]

Reimplemented from SyFi::StandardFE.

Definition at line 2591 of file SyFi.py.

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

Reimplemented from SyFi::StandardFE.

Definition at line 2575 of file SyFi.py.

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

Reimplemented from SyFi::StandardFE.

Definition at line 2572 of file SyFi.py.

Definition at line 15 of file Nedelec2Hdiv.h.

Referenced by compute_basis_functions().

Reimplemented from SyFi::StandardFE.

Definition at line 2590 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