|
SyFi
0.3
|
|
SyFi
0.3
|
#include <BrezziDouglasMarini.h>
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 |
Definition at line 26 of file BrezziDouglasMarini.h.
Definition at line 30 of file BrezziDouglasMarini.cpp.
References SyFi::StandardFE::description.
: StandardFE() { description = "BrezziDouglasMarini"; }
| SyFi::BrezziDouglasMarini::BrezziDouglasMarini | ( | Polygon & | p, |
| int | order = 1, |
||
| bool | pointwise = true |
||
| ) |
Definition at line 35 of file BrezziDouglasMarini.cpp.
References compute_basis_functions(), and pointwise.
: StandardFE(p, order) { pointwise = pointwise_; compute_basis_functions(); }
| virtual SyFi::BrezziDouglasMarini::~BrezziDouglasMarini | ( | ) | [inline, virtual] |
Definition at line 33 of file BrezziDouglasMarini.h.
{}
| void SyFi::BrezziDouglasMarini::compute_basis_functions | ( | ) | [virtual] |
Reimplemented from SyFi::StandardFE.
Definition at line 41 of file BrezziDouglasMarini.cpp.
References SyFi::bernsteinv(), test_syfi::debug::c, 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, SyFi::x, and SyFi::y.
Referenced by BrezziDouglasMarini().
{
if ( order < 1 )
{
throw(std::logic_error("Brezzi-Douglas-Marini elements must be of order 1 or higher."));
}
if ( p == NULL )
{
throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
}
GiNaC::ex nsymb = GiNaC::symbol("n");
if (pointwise)
{
if ( p->str().find("ReferenceLine") != string::npos )
{
cout <<"Can not define the Brezzi-Douglas-Marini element on a line"<<endl;
}
if ( p->str().find("ReferenceTetrahedron") != string::npos )
{
cout <<"Can not define the Brezzi-Douglas-Marini element on a tetrahedron, see Nedelec2Hdiv"<<endl;
}
else if ( p->str().find("Triangle") != string::npos )
{
description = istr("BrezziDouglasMarini_", order) + "_2D";
Triangle& triangle = (Triangle&)(*p);
GiNaC::lst equations;
GiNaC::lst variables;
GiNaC::lst polynom_space = bernsteinv(2,order, triangle, "b");
GiNaC::ex pspace = polynom_space.op(0);
variables = collapse(GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)));
GiNaC::symbol t("t");
GiNaC::ex dofi;
// dofs related to edges
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::lst normal_vec = normal(triangle, i);
GiNaC::ex Vn = inner(pspace, normal_vec);
GiNaC::lst points = interior_coordinates(line, order);
GiNaC::ex point;
for (unsigned int j=0; j< points.nops(); j++)
{
point = points.op(j);
dofi = Vn.subs(x == point.op(0)).subs(y == point.op(1));
GiNaC::ex eq = dofi == GiNaC::numeric(0);
equations.append(eq);
dofs.insert(dofs.end(), GiNaC::lst(point,nsymb));
}
}
// std::cout <<"no variables "<<variables.nops()<<std::endl;
// std::cout <<"no equations "<<equations.nops()<<std::endl;
// invert the matrix:
// GiNaC has a bit strange way to invert a matrix.
// It solves the system AA^{-1} = Id.
// It seems that this way is the only way to do
// properly with the solve_algo::gauss flag.
//
GiNaC::matrix b; GiNaC::matrix A;
matrix_from_equations(equations, variables, A, b);
unsigned int ncols = A.cols();
GiNaC::matrix vars_sq(ncols, ncols);
// matrix of symbols
for (unsigned r=0; r<ncols; ++r)
for (unsigned c=0; c<ncols; ++c)
vars_sq(r, c) = GiNaC::symbol();
GiNaC::matrix id(ncols, ncols);
// identity
const GiNaC::ex _ex1(1);
for (unsigned i=0; i<ncols; ++i)
id(i, i) = _ex1;
// invert the matrix
GiNaC::matrix m_inv(ncols, ncols);
m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss);
for (unsigned int i=0; i<dofs.size(); i++)
{
b.let_op(i) = GiNaC::numeric(1);
GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b));
GiNaC::lst subs;
for (unsigned int ii=0; ii<xx.nops(); ii++)
{
subs.append(variables.op(ii) == xx.op(ii));
}
GiNaC::ex Nj1 = pspace.op(0).subs(subs);
GiNaC::ex Nj2 = pspace.op(1).subs(subs);
Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2)));
b.let_op(i) = GiNaC::numeric(0);
}
}
}
}
| GiNaC::lst SyFi::BrezziDouglasMarini::dof_repr |
Definition at line 30 of file BrezziDouglasMarini.h.
Definition at line 29 of file BrezziDouglasMarini.h.
Referenced by BrezziDouglasMarini(), and compute_basis_functions().
1.7.6.1