SyFi::BrezziDouglasMarini Class Reference

#include <BrezziDouglasMarini.h>

Inheritance diagram for SyFi::BrezziDouglasMarini:

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

List of all members.

Public Member Functions

 BrezziDouglasMarini ()
 BrezziDouglasMarini (Polygon &p, int order=1, bool pointwise=true)
virtual ~BrezziDouglasMarini ()
virtual void compute_basis_functions ()

Public Attributes

bool pointwise
GiNaC::lst dof_repr


Detailed Description

Definition at line 12 of file BrezziDouglasMarini.h.


Constructor & Destructor Documentation

SyFi::BrezziDouglasMarini::BrezziDouglasMarini (  ) 

Definition at line 16 of file BrezziDouglasMarini.cpp.

References SyFi::StandardFE::description.

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

SyFi::BrezziDouglasMarini::BrezziDouglasMarini ( Polygon p,
int  order = 1,
bool  pointwise = true 
)

Definition at line 21 of file BrezziDouglasMarini.cpp.

References compute_basis_functions(), and pointwise.

00021                                                                                         : StandardFE(p, order)
00022         {
00023                 pointwise = pointwise_;
00024                 compute_basis_functions();
00025         }

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

Definition at line 19 of file BrezziDouglasMarini.h.

00019 {}


Member Function Documentation

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

Reimplemented from SyFi::StandardFE.

Definition at line 27 of file BrezziDouglasMarini.cpp.

References SyFi::bernsteinv(), SyFi::collapse(), SyFi::StandardFE::description, SyFi::StandardFE::dofs, SyFi::inner(), SyFi::interior_coordinates(), SyFi::istr(), SyFi::Triangle::line(), SyFi::matrix_from_equations(), SyFi::normal(), SyFi::StandardFE::Ns, SyFi::StandardFE::order, SyFi::StandardFE::p, pointwise, SyFi::Polygon::str(), SyFi::t, demos::crouzeixraviart::triangle, SyFi::x, and SyFi::y.

Referenced by BrezziDouglasMarini().

00028         {
00029 
00030                 if ( order < 1 )
00031                 {
00032                         throw(std::logic_error("Brezzi-Douglas-Marini elements must be of order 1 or higher."));
00033                 }
00034 
00035                 if ( p == NULL )
00036                 {
00037                         throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00038                 }
00039 
00040                 GiNaC::ex nsymb = GiNaC::symbol("n");
00041                 if (pointwise)
00042                 {
00043 
00044                         if ( p->str().find("ReferenceLine") != string::npos )
00045                         {
00046 
00047                                 cout <<"Can not define the Brezzi-Douglas-Marini element on a line"<<endl;
00048 
00049                         }
00050 
00051                         if ( p->str().find("ReferenceTetrahedron") != string::npos )
00052                         {
00053 
00054                                 cout <<"Can not define the Brezzi-Douglas-Marini element on a tetrahedron, see Nedelec2Hdiv"<<endl;
00055 
00056                         }
00057                         else if ( p->str().find("Triangle") != string::npos )
00058                         {
00059 
00060                                 description = istr("BrezziDouglasMarini_", order) + "_2D";
00061 
00062                                 Triangle& triangle = (Triangle&)(*p);
00063                                 GiNaC::lst equations;
00064                                 GiNaC::lst variables;
00065                                 GiNaC::lst polynom_space = bernsteinv(2,order, triangle, "b");
00066                                 GiNaC::ex pspace = polynom_space.op(0);
00067 
00068                                 variables = collapse(GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)));
00069 
00070                                 GiNaC::symbol t("t");
00071                                 GiNaC::ex dofi;
00072                                 // dofs related to edges
00073                                 for (int i=0; i< 3; i++)
00074                                 {
00075                                         Line line = triangle.line(i);
00076                                         GiNaC::lst normal_vec = normal(triangle, i);
00077                                         GiNaC::ex Vn = inner(pspace, normal_vec);
00078                                         GiNaC::lst points = interior_coordinates(line, order);
00079 
00080                                         GiNaC::ex point;
00081                                         for (unsigned int j=0; j< points.nops(); j++)
00082                                         {
00083                                                 point = points.op(j);
00084                                                 dofi = Vn.subs(x == point.op(0)).subs(y == point.op(1));
00085                                                 GiNaC::ex eq = dofi == GiNaC::numeric(0);
00086                                                 equations.append(eq);
00087                                                 dofs.insert(dofs.end(), GiNaC::lst(point,nsymb));
00088                                         }
00089                                 }
00090 
00091 //                              std::cout <<"no variables "<<variables.nops()<<std::endl;
00092 //                              std::cout <<"no equations "<<equations.nops()<<std::endl;
00093 
00094                                 // invert the matrix:
00095                                 // GiNaC has a bit strange way to invert a matrix.
00096                                 // It solves the system AA^{-1} = Id.
00097                                 // It seems that this way is the only way to do
00098                                 // properly with the solve_algo::gauss flag.
00099                                 //
00100                                 GiNaC::matrix b; GiNaC::matrix A;
00101                                 matrix_from_equations(equations, variables, A, b);
00102 
00103                                 unsigned int ncols = A.cols();
00104                                 GiNaC::matrix vars_sq(ncols, ncols);
00105 
00106                                 // matrix of symbols
00107                                 for (unsigned r=0; r<ncols; ++r)
00108                                         for (unsigned c=0; c<ncols; ++c)
00109                                                 vars_sq(r, c) = GiNaC::symbol();
00110 
00111                                 GiNaC::matrix id(ncols, ncols);
00112 
00113                                 // identity
00114                                 const GiNaC::ex _ex1(1);
00115                                 for (unsigned i=0; i<ncols; ++i)
00116                                         id(i, i) = _ex1;
00117 
00118                                 // invert the matrix
00119                                 GiNaC::matrix m_inv(ncols, ncols);
00120                                 m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss);
00121 
00122                                 for (unsigned int i=0; i<dofs.size(); i++)
00123                                 {
00124                                         b.let_op(i) = GiNaC::numeric(1);
00125                                         GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b));
00126 
00127                                         GiNaC::lst subs;
00128                                         for (unsigned int ii=0; ii<xx.nops(); ii++)
00129                                         {
00130                                                 subs.append(variables.op(ii) == xx.op(ii));
00131                                         }
00132                                         GiNaC::ex Nj1 = pspace.op(0).subs(subs);
00133                                         GiNaC::ex Nj2 = pspace.op(1).subs(subs);
00134                                         Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2)));
00135                                         b.let_op(i) = GiNaC::numeric(0);
00136                                 }
00137                         }
00138                 }
00139 
00140         }


Member Data Documentation

Definition at line 16 of file BrezziDouglasMarini.h.

Definition at line 15 of file BrezziDouglasMarini.h.

Referenced by BrezziDouglasMarini(), and compute_basis_functions().


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

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