|
SyFi
0.3
|
|
SyFi
0.3
|
#include <CrouzeixRaviart.h>
Public Member Functions | |
| CrouzeixRaviart () | |
| CrouzeixRaviart (Polygon &p, unsigned int order=1) | |
| virtual | ~CrouzeixRaviart () |
| void | compute_basis_functions () |
Definition at line 26 of file CrouzeixRaviart.h.
Definition at line 29 of file CrouzeixRaviart.cpp.
References SyFi::StandardFE::order.
: StandardFE() { order = 1; }
| SyFi::CrouzeixRaviart::CrouzeixRaviart | ( | Polygon & | p, |
| unsigned int | order = 1 |
||
| ) |
Definition at line 34 of file CrouzeixRaviart.cpp.
References compute_basis_functions().
: StandardFE(p, order) { compute_basis_functions(); }
| virtual SyFi::CrouzeixRaviart::~CrouzeixRaviart | ( | ) | [inline, virtual] |
Definition at line 31 of file CrouzeixRaviart.h.
{}
| void SyFi::CrouzeixRaviart::compute_basis_functions | ( | ) | [virtual] |
Reimplemented from SyFi::StandardFE.
Definition at line 39 of file CrouzeixRaviart.cpp.
References SyFi::bernstein(), SyFi::StandardFE::description, SyFi::dirac(), SyFi::StandardFE::dofs, SyFi::Line::integrate(), SyFi::Triangle::line(), SyFi::StandardFE::Ns, SyFi::StandardFE::order, SyFi::StandardFE::p, SyFi::Polygon::str(), SyFi::sub(), SyFi::t, SyFi::Tetrahedron::triangle(), SyFi::Polygon::vertex(), SyFi::x, SyFi::y, and SyFi::z.
Referenced by check_CrouzeixRaviart(), SyFi::VectorCrouzeixRaviart::compute_basis_functions(), and CrouzeixRaviart().
{
// remove previously computed basis functions and dofs
Ns.clear();
dofs.clear();
if ( p == NULL )
{
throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
}
if (order != 1)
{
throw(std::logic_error("Only Crouziex-Raviart elements of order 1 is possible"));
}
// see e.g. Brezzi and Fortin book page 116 for the definition
if ( p->str().find("ReferenceLine") != string::npos )
{
cout <<"Can not define the Raviart-Thomas element on a line"<<endl;
}
else if ( p->str().find("Triangle") != string::npos )
{
description = "CrouzeixRaviart_2D";
Triangle& triangle = (Triangle&)(*p);
// create the polynomial space
GiNaC::ex polynom_space = bernstein(1, triangle, "a");
GiNaC::ex polynom = polynom_space.op(0);
GiNaC::lst variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
GiNaC::ex basis = polynom_space.op(2);
// create the dofs
GiNaC::symbol t("t");
for (int j=0; j< 3; j++)
{
// solve the linear system to compute
// each of the basis functions
GiNaC::lst equations;
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
// GiNaC::ex dofi = line.integrate(polynom);
GiNaC::lst midpoint = GiNaC::lst(
(line.vertex(0).op(0) + line.vertex(1).op(0))/2,
(line.vertex(0).op(1) + line.vertex(1).op(1))/2);
dofs.insert(dofs.end(), midpoint);
// GiNaC::ex dofi = polynom.subs( x == midpoint.op(0)).subs( y == midpoint.op(1));
GiNaC::ex dofi = line.integrate(polynom);
equations.append( dofi == dirac(i,j));
if (j == 1)
{
// GiNaC::lst d = GiNaC::lst(line.vertex(0) , line.vertex(1));
dofs.insert(dofs.end(), midpoint);
}
}
GiNaC::ex sub = lsolve(equations, variables);
GiNaC::ex Ni = polynom.subs(sub);
Ns.insert(Ns.end(),Ni);
}
}
else if ( p->str().find("Tetrahedron") != string::npos )
{
description = "CrouzeixRaviart_3D";
Tetrahedron& tetrahedron = (Tetrahedron&)(*p);
GiNaC::ex polynom_space = bernstein(1, tetrahedron, "a");
GiNaC::ex polynom = polynom_space.op(0);
GiNaC::lst variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
GiNaC::ex basis = polynom_space.op(2);
GiNaC::ex bernstein_pol;
GiNaC::symbol t("t");
// dofs related to edges
for (int j=0; j< 4; j++)
{
GiNaC::lst equations;
for (int i=0; i< 4; i++)
{
Triangle triangle = tetrahedron.triangle(i);
GiNaC::lst midpoint = GiNaC::lst(
(triangle.vertex(0).op(0) + triangle.vertex(1).op(0) + triangle.vertex(2).op(0))/3,
(triangle.vertex(0).op(1) + triangle.vertex(1).op(1) + triangle.vertex(2).op(1))/3,
(triangle.vertex(0).op(2) + triangle.vertex(1).op(2) + triangle.vertex(2).op(2))/3
);
GiNaC::ex dofi = polynom.subs(x == midpoint.op(0)).subs(y == midpoint.op(1)).subs(z == midpoint.op(2));
equations.append( dofi == dirac(i,j));
if ( j == 1 )
{
// GiNaC::lst d = GiNaC::lst(triangle.vertex(0), triangle.vertex(1), triangle.vertex(2));
dofs.insert(dofs.end(), midpoint);
}
}
GiNaC::ex sub = lsolve(equations, variables);
GiNaC::ex Ni = polynom.subs(sub);
Ns.insert(Ns.end(),Ni);
}
}
}
1.7.6.1