|
SyFi
0.3
|
|
SyFi
0.3
|
#include <Lagrange.h>
Public Member Functions | |
| Lagrange () | |
| Lagrange (Polygon &p, unsigned int order=1) | |
| virtual | ~Lagrange () |
| virtual void | compute_basis_functions () |
Definition at line 26 of file Lagrange.h.
Definition at line 31 of file Lagrange.cpp.
References SyFi::StandardFE::description.
: StandardFE() { description = "Lagrange"; }
| SyFi::Lagrange::Lagrange | ( | Polygon & | p, |
| unsigned int | order = 1 |
||
| ) |
Definition at line 36 of file Lagrange.cpp.
References compute_basis_functions().
: StandardFE(p, order) { compute_basis_functions(); }
| virtual SyFi::Lagrange::~Lagrange | ( | ) | [inline, virtual] |
Definition at line 31 of file Lagrange.h.
{}
| void SyFi::Lagrange::compute_basis_functions | ( | ) | [virtual] |
Reimplemented from SyFi::StandardFE.
Reimplemented in SyFi::DiscontinuousLagrange.
Definition at line 41 of file Lagrange.cpp.
References SyFi::bernstein(), SyFi::bezier_ordinates(), test_syfi::debug::c, SyFi::StandardFE::description, SyFi::StandardFE::dof(), SyFi::StandardFE::dofs, SyFi::istr(), SyFi::matrix_from_equations(), SyFi::StandardFE::N(), SyFi::StandardFE::nbf(), SyFi::StandardFE::Ns, SyFi::StandardFE::order, SyFi::StandardFE::p, SyFi::pol(), SyFi::Polygon::str(), SyFi::t, SyFi::x, SyFi::y, and SyFi::z.
Referenced by SyFi::VectorLagrange::compute_basis_functions(), SyFi::TensorLagrange::compute_basis_functions(), Lagrange(), and main().
{
// NOTE: in the below code dof(i) is not used to
// determine the basis functions
// remove previously computed basis functions and dofs
Ns.clear();
dofs.clear();
if ( order < 1 )
{
throw(std::logic_error("Lagrangian 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::lst equations;
GiNaC::lst variables;
GiNaC::ex polynom;
if (p->str().find("ReferenceLine") != string::npos)
{
description = istr("Lagrange_", order) + "_1D";
// Look at the case with the Triangle for a documented code
// polynom = pol(order, 1, "a");
// variables = coeffs(polynom);
GiNaC::ex polynom_space = bernstein(order, *p, "a");
// GiNaC::ex polynom_space = pol(order, 1, "a");
polynom = polynom_space.op(0);
variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
GiNaC::ex increment = GiNaC::numeric(1,order);
GiNaC::ex Nj;
for (GiNaC::ex p=0; p<= 1 ; p += increment )
{
GiNaC::ex eq = polynom == GiNaC::numeric(0);
equations.append(eq.subs(x == p));
dofs.insert(dofs.end(), GiNaC::lst(p));
}
}
else if (p->str().find("Line") != string::npos )
{
description = istr("Lagrange_", order) + "_1D";
// Look at the case with the Triangle for a documented code
// polynom = pol(order, 1, "a");
// variables = coeffs(polynom);
GiNaC::ex polynom_space = bernstein(order, *p, "a");
// GiNaC::ex polynom_space = pol(order, 1, "a");
polynom = polynom_space.op(0);
variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
Polygon& pp = *p;
Line& l = (Line&) pp;
GiNaC::lst points = bezier_ordinates(l,order);
GiNaC::ex Nj;
for (unsigned int i=1; i<= points.nops() ; i++ )
{
GiNaC::ex point = points.op(i-1);
GiNaC::ex eq = polynom == GiNaC::numeric(0);
if (point.nops() == 0) eq = eq.subs(x == point);
if (point.nops() > 0) eq = eq.subs(x == point.op(0));
if (point.nops() > 1) eq = eq.subs(y == point.op(1));
if (point.nops() > 2) eq = eq.subs(z == point.op(2));
equations.append(eq);
dofs.insert(dofs.end(), GiNaC::lst(point));
}
}
else if (p->str().find("ReferenceTriangle") != string::npos )
{
description = istr("Lagrange_", order) + "_2D";
// Look at the case with the Triangle for a documented code
// polynom = pol(order, 2, "b");
// variables = coeffs(polynom);
GiNaC::ex polynom_space = bernstein(order, *p, "b");
// GiNaC::ex polynom_space = pol(order, 2, "a");
polynom = polynom_space.op(0);
variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
GiNaC::ex increment = GiNaC::numeric(1,order);
GiNaC::ex Nj;
GiNaC::numeric one = 1;
for (GiNaC::ex q=0; q<= one ; q += increment )
{
for (GiNaC::ex p=0; p<= one-q ; p += increment )
{
GiNaC::ex eq = polynom == GiNaC::numeric(0);
equations.append(eq.subs(GiNaC::lst(x == p, y == q)));
dofs.insert(dofs.end(), GiNaC::lst(p,q));
}
}
}
else if ( p->str().find("Triangle") != string::npos)
{
description = istr("Lagrange_", order) + "_2D";
// Look HERE for the documented code
// GiNaC::ex polynom_space = pol(order, 2, "a");
GiNaC::ex polynom_space = bernstein(order, *p, "b");
// the polynomial spaces on the form:
// first item: a0 + a1*x + a2*y + a3*x^2 + a4*x*y ... the polynom
// second item: a0, a1, a2, ... the coefficents
// third item 1, x, y, x^2, .. the basis
// Could also do:
// GiNaC::ex polynom_space = bernstein(order, t, "a");
polynom = polynom_space.op(0);
variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
GiNaC::ex Nj;
Polygon& pp = *p;
Triangle& t = (Triangle&) pp;
// The bezier ordinates (in which the basis function should be either 0 or 1)
GiNaC::lst points = bezier_ordinates(t,order);
for (unsigned int i=1; i<= points.nops() ; i++ )
{
GiNaC::ex point = points.op(i-1);
GiNaC::ex eq = polynom == GiNaC::numeric(0);
equations.append(eq.subs(GiNaC::lst(x == point.op(0) , y == point.op(1))));
dofs.insert(dofs.end(), GiNaC::lst(point.op(0),point.op(1)));
}
}
else if ( p->str().find("ReferenceRectangle") != string::npos)
{
description = istr("Lagrange_", order) + "_2D";
// create 1D element, then create tensor product
ReferenceLine line;
Lagrange fe(line, order);
for (unsigned int i=0; i< fe.nbf(); i++)
{
for (unsigned int j=0; j< fe.nbf(); j++)
{
Ns.insert(Ns.end(), fe.N(i)*fe.N(j).subs(x==y));
dofs.insert(dofs.end(), GiNaC::lst(fe.dof(i).op(0), fe.dof(j).op(0)));
}
}
return;
/* OLD CODE
GiNaC::ex polynom_space = legendre(order, 2, "b");
polynom = polynom_space.op(0);
variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
GiNaC::ex Nj;
Polygon& pp = *p;
ReferenceRectangle& b = (ReferenceRectangle&) pp;
int no_points = (order+1)*(order+1);
GiNaC::ex increment = GiNaC::numeric(1,order);
GiNaC::numeric one=1.0;
for (GiNaC::ex q=0; q <= one ; q += increment )
{
for (GiNaC::ex p=0; p <= one ; p += increment )
{
GiNaC::ex eq = polynom == GiNaC::numeric(0);
equations.append(eq.subs(GiNaC::lst(x == p, y == q)));
dofs.push_back(GiNaC::lst(p,q));
}
}
*/
}
else if ( p->str().find("ReferenceTetrahedron") != string::npos)
{
description = istr("Lagrange_", order) + "_3D";
// Look at the case with the Triangle for a documented code
// polynom = pol(order, 3, "b");
// GiNaC::ex polynom_space = pol(order, 3, "a");
GiNaC::ex polynom_space = bernstein(order, *p, "b");
polynom = polynom_space.op(0);
variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
int nno =0;
for (unsigned int j=0; j<= order; j++)
{
nno += (j+1)*(j+2)/2;
}
GiNaC::ex increment = GiNaC::numeric(1,order);
GiNaC::ex Nj;
for (GiNaC::ex r=0; r<= 1 ; r += increment )
{
for (GiNaC::ex q=0; q<= 1-r ; q += increment )
{
for (GiNaC::ex p=0; p<= 1-r-q ; p += increment )
{
GiNaC::ex eq = polynom == GiNaC::numeric(0);
equations.append(eq.subs(GiNaC::lst(x == p, y == q, z == r )));
dofs.push_back(GiNaC::lst(p,q,r));
}
}
}
}
else if ( p->str().find("Tetrahedron") != string::npos)
{
description = istr("Lagrange_", order) + "_3D";
// Look at the case with the Triangle for a documented code
GiNaC::ex polynom_space = pol(order, 3, "a");
// GiNaC::ex polynom_space = bernstein(order, *p, "b");
polynom = polynom_space.op(0);
variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
GiNaC::ex increment = GiNaC::numeric(1,order);
GiNaC::ex Nj;
Polygon& pp = *p;
Tetrahedron& t = (Tetrahedron&) pp;
GiNaC::lst points = bezier_ordinates(t,order);
for (unsigned int i=1; i<= points.nops() ; i++ )
{
GiNaC::ex point = points.op(i-1);
GiNaC::ex eq = polynom == GiNaC::numeric(0);
equations.append(eq.subs(GiNaC::lst(x == point.op(0) , y == point.op(1), z == point.op(2))));
dofs.push_back(GiNaC::lst(point.op(0),point.op(1),point.op(2)));
}
}
else if ( p->str().find("ReferenceBox") != string::npos)
{
description = istr("Lagrange_", order) + "_3D";
ReferenceLine line;
Lagrange fe(line, order);
for (unsigned int i=0; i< fe.nbf(); i++)
{
for (unsigned int j=0; j< fe.nbf(); j++)
{
for (unsigned int k=0; k< fe.nbf(); k++)
{
Ns.insert(Ns.end(), fe.N(i)*fe.N(j).subs(x==y)*fe.N(k).subs(x==z));
dofs.insert(dofs.end(), GiNaC::lst(fe.dof(i).op(0), fe.dof(j).op(0), fe.dof(k).op(0)));
}
}
}
return;
/* OLD CODE
GiNaC::ex polynom_space = legendre(order, 3, "b");
polynom = polynom_space.op(0);
variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
GiNaC::ex Nj;
Polygon& pp = *p;
ReferenceRectangle& b = (ReferenceRectangle&) pp;
int no_points = (order+1)*(order+1)*(order+1);
GiNaC::ex increment = GiNaC::numeric(1,order);
GiNaC::numeric one = 1;
for (GiNaC::ex p=0; p <= one ; p += increment) {
for (GiNaC::ex q=0; q <= one ; q += increment) {
for (GiNaC::ex r=0; r <= one ; r += increment) {
GiNaC::ex eq = polynom == GiNaC::numeric(0);
equations.append(eq.subs(GiNaC::lst(x == p, y == q, z==r)));
dofs.push_back(GiNaC::lst(p,q,r));
}
}
}
/ *
GiNaC::ex subs = lsolve(equations, variables);
Nj = polynom.subs(subs);
Ns.push_back(Nj);
*/
}
// 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.
//
// std::cout <<"no variables "<<variables.nops()<<std::endl;
// std::cout <<"no equations "<<equations.nops()<<std::endl;
// print(equations);
// print(variables);
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 Nj = polynom.subs(subs);
Ns.insert(Ns.end(), Nj);
b.let_op(i) = GiNaC::numeric(0);
}
}
1.7.6.1