|
SyFi
0.3
|
|
SyFi
0.3
|
#include <Robust.h>
Public Member Functions | |
| Robust () | |
| Robust (Polygon &p, unsigned int order=0, bool pointwise=true) | |
| virtual | ~Robust () |
| virtual void | compute_basis_functions () |
| virtual void | compute_basis_functions_old () |
Public Attributes | |
| bool | pointwise |
| GiNaC::lst | dof_repr |
Definition at line 30 of file Robust.cpp.
References SyFi::StandardFE::description.
: StandardFE() { description = "Robust"; }
| SyFi::Robust::Robust | ( | Polygon & | p, |
| unsigned int | order = 0, |
||
| bool | pointwise = true |
||
| ) |
| virtual SyFi::Robust::~Robust | ( | ) | [inline, virtual] |
| void SyFi::Robust::compute_basis_functions | ( | ) | [virtual] |
Reimplemented from SyFi::StandardFE.
Definition at line 41 of file Robust.cpp.
References SyFi::bernstein(), SyFi::bernsteinv(), test_syfi::debug::c, SyFi::coeff(), SyFi::collapse(), SyFi::StandardFE::description, SyFi::div(), dof_repr, SyFi::StandardFE::dofs, SyFi::inner(), SyFi::Line::integrate(), SyFi::Triangle::integrate(), 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::pol2basisandcoeff(), SyFi::Line::repr(), SyFi::Polygon::str(), SyFi::t, SyFi::tangent(), SyFi::Polygon::vertex(), SyFi::x, and SyFi::y.
Referenced by main().
{
// remove previously computed basis functions and dofs
Ns.clear();
dofs.clear();
GiNaC::ex nsymb = GiNaC::symbol("n");
GiNaC::ex tsymb = GiNaC::symbol("t");
if ( p == NULL )
{
throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
}
if ( p->str().find("Line") != string::npos )
{
cout <<"Can not define the Robust element on a line"<<endl;
}
else if ( p->str().find("Triangle") != string::npos )
{
if (pointwise)
{
description = "Robust_2D";
Triangle& triangle = (Triangle&)(*p);
GiNaC::lst equations;
GiNaC::lst variables;
GiNaC::ex V_space = bernsteinv(2, order+3, triangle, "a");
GiNaC::ex V = V_space.op(0);
GiNaC::ex V_vars = V_space.op(1);
GiNaC::ex divV = div(V);
exmap basis2coeff = pol2basisandcoeff(divV);
exmap::iterator iter;
// div constraints:
for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
{
GiNaC::ex basis = (*iter).first;
GiNaC::ex coeff= (*iter).second;
if ( coeff != 0 )
{
}
if ( coeff != 0 && ( basis.degree(x) + basis.degree(y) > int(order) ) )
{
equations.append( coeff == 0 );
}
}
// normal constraints on edges:
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::symbol s("s");
GiNaC::lst normal_vec = normal(triangle, i);
GiNaC::ex Vn = inner(V, normal_vec);
Vn = Vn.subs(line.repr(s).op(0)).subs(line.repr(s).op(1));
basis2coeff = pol2basisandcoeff(Vn,s);
for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
{
GiNaC::ex basis = (*iter).first;
GiNaC::ex coeff= (*iter).second;
if ( coeff != 0 && basis.degree(s) > int(order+1) )
{
equations.append( coeff == 0 );
}
}
}
if ( order%2==1 )
{
// tangent constraints on edges:
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::symbol s("s");
GiNaC::lst tangent_vec = tangent(triangle, i);
GiNaC::ex Vt = inner(V, tangent_vec);
Vt = Vt.subs(line.repr(s).op(0)).subs(line.repr(s).op(1));
basis2coeff = pol2basisandcoeff(Vt,s);
for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
{
GiNaC::ex basis = (*iter).first;
GiNaC::ex coeff= (*iter).second;
if ( coeff != 0 && basis.degree(s) > int(order+2) )
{
equations.append( coeff == 0 );
}
}
}
}
GiNaC::ex dofi;
// dofs related to the normal on the edges
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::lst normal_vec = normal(triangle, i);
GiNaC::ex Vn = inner(V, normal_vec);
GiNaC::lst points = interior_coordinates(line, order + 1);
GiNaC::ex edge_length = line.integrate(GiNaC::numeric(1));
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))*edge_length;
dofs.insert(dofs.end(), GiNaC::lst(point,nsymb));
GiNaC::ex eq = dofi == GiNaC::numeric(0);
equations.append(eq);
}
}
if ( order%2==0)
{
// dofs related to the tangent on the edges
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::lst tangent_vec = tangent(triangle, i);
GiNaC::ex Vt = inner(V, tangent_vec);
GiNaC::lst points = interior_coordinates(line, order);
GiNaC::ex edge_length = line.integrate(GiNaC::numeric(1));
GiNaC::ex point;
for (unsigned int j=0; j< points.nops(); j++)
{
point = points.op(j);
dofi = Vt.subs(x == point.op(0)).subs(y == point.op(1))*edge_length;
dofs.insert(dofs.end(), GiNaC::lst(point,tsymb));
GiNaC::ex eq = dofi == GiNaC::numeric(0);
equations.append(eq);
}
}
}
if ( order%2==1 )
{
// dofs related to the tangent on the edges
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::lst tangent_vec = tangent(triangle, i);
GiNaC::ex Vt = inner(V, tangent_vec);
GiNaC::lst points = interior_coordinates(line, order-1);
GiNaC::ex edge_length = line.integrate(GiNaC::numeric(1));
GiNaC::ex point;
for (unsigned int j=0; j< points.nops(); j++)
{
point = points.op(j);
dofi = Vt.subs(x == point.op(0)).subs(y == point.op(1))*edge_length;
dofs.insert(dofs.end(), GiNaC::lst(point,tsymb));
GiNaC::ex eq = dofi == GiNaC::numeric(0);
equations.append(eq);
}
}
}
// dofs related to the whole triangle
GiNaC::lst bernstein_polv;
if ( order > 0)
{
GiNaC::lst points = interior_coordinates(triangle, order-1);
GiNaC::ex point;
GiNaC::ex eq;
for (unsigned int i=0; i< points.nops(); i++)
{
point = points.op(i);
// x - components of interior dofs
dofi = V.op(0).subs(x == point.op(0)).subs(y == point.op(1));
eq = dofi == GiNaC::numeric(0);
equations.append(eq);
dofs.insert(dofs.end(), GiNaC::lst(point,0));
// y - components of interior dofs
dofi = V.op(1).subs(x == point.op(0)).subs(y == point.op(1));
eq = dofi == GiNaC::numeric(0);
equations.append(eq);
dofs.insert(dofs.end(), GiNaC::lst(point,1));
}
}
variables = collapse(GiNaC::lst(V_vars.op(0),V_vars.op(1)));
// std::cout <<"no variables "<<variables.nops()<<std::endl;
// std::cout <<"no equations "<<equations.nops()<<std::endl;
GiNaC::matrix b; GiNaC::matrix A;
matrix_from_equations(equations, variables, A, b);
// std::cout <<" A "<<A.evalf()<<std::endl;
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(11+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 = V.op(0).subs(subs);
GiNaC::ex Nj2 = V.op(1).subs(subs);
Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2)));
b.let_op(11+i) = GiNaC::numeric(0);
}
} //dofs are integrals etc. ie not represented as normal/tangents in points
else
{
description = "Robust_2D";
Triangle& triangle = (Triangle&)(*p);
GiNaC::lst equations;
GiNaC::lst variables;
GiNaC::ex V_space = bernsteinv(2, order+3, triangle, "a");
GiNaC::ex V = V_space.op(0);
GiNaC::ex V_vars = V_space.op(1);
GiNaC::ex divV = div(V);
exmap basis2coeff = pol2basisandcoeff(divV);
exmap::iterator iter;
// div constraints:
for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
{
GiNaC::ex basis = (*iter).first;
GiNaC::ex coeff= (*iter).second;
if ( coeff != 0 )
{
}
if ( coeff != 0 && ( basis.degree(x) + basis.degree(y) > int(order) ) )
{
equations.append( coeff == 0 );
}
}
// normal constraints on edges:
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::symbol s("s");
GiNaC::lst normal_vec = normal(triangle, i);
GiNaC::ex Vn = inner(V, normal_vec);
Vn = Vn.subs(line.repr(s).op(0)).subs(line.repr(s).op(1));
basis2coeff = pol2basisandcoeff(Vn,s);
for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
{
GiNaC::ex basis = (*iter).first;
GiNaC::ex coeff= (*iter).second;
if ( coeff != 0 && basis.degree(s) > int(order+1) )
{
equations.append( coeff == 0 );
}
}
}
if ( order%2==1 )
{
// tangent constraints on edges:
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::symbol s("s");
GiNaC::lst tangent_vec = tangent(triangle, i);
GiNaC::ex Vt = inner(V, tangent_vec);
Vt = Vt.subs(line.repr(s).op(0)).subs(line.repr(s).op(1));
basis2coeff = pol2basisandcoeff(Vt,s);
for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
{
GiNaC::ex basis = (*iter).first;
GiNaC::ex coeff= (*iter).second;
if ( coeff != 0 && basis.degree(s) > int(order+2) )
{
equations.append( coeff == 0 );
}
}
}
}
// dofs related to the normal on the edges
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::lst normal_vec = normal(triangle, i);
GiNaC::ex Pk1_space = bernstein(order+1, line, istr("a",i));
GiNaC::ex Pk1 = Pk1_space.op(2);
GiNaC::ex Vn = inner(V, normal_vec);
GiNaC::ex basis;
for (unsigned int j=0; j< Pk1.nops(); j++)
{
basis = Pk1.op(j);
GiNaC::ex integrand = Vn*basis;
GiNaC::ex dofi = line.integrate(integrand);
// dofs.insert(dofs.end(), dofi);
GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),j);
dofs.insert(dofs.end(), d);
GiNaC::ex eq = dofi == GiNaC::numeric(0);
equations.append(eq);
GiNaC::ex u = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]")));
GiNaC::ex n = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("normal_vec[0]"), GiNaC::symbol("normal_vec[1]")));
dof_repr.append(GiNaC::lst(inner(u,n)*basis.subs( x == GiNaC::symbol("xi[0]"))
.subs( y == GiNaC::symbol("xi[1]")), d));
}
}
// dofs related to the tangent on the edges
if ( order%2==0 )
{
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::lst tangent_vec = tangent(triangle, i);
GiNaC::ex Pk_space = bernstein(order, line, istr("a",i));
GiNaC::ex Pk = Pk_space.op(2);
GiNaC::ex Vt = inner(V, tangent_vec);
GiNaC::ex basis;
for (unsigned int j=0; j< Pk.nops(); j++)
{
basis = Pk.op(j);
GiNaC::ex integrand = Vt*basis;
GiNaC::ex dofi = line.integrate(integrand);
// dofs.insert(dofs.end(), dofi);
GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),2);
dofs.insert(dofs.end(), d);
GiNaC::ex eq = dofi == GiNaC::numeric(0);
equations.append(eq);
GiNaC::ex u = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]")));
GiNaC::ex t = GiNaC::matrix(2,1,GiNaC::lst(-GiNaC::symbol("normal_vec[1]"), GiNaC::symbol("normal_vec[0]")));
dof_repr.append(GiNaC::lst(inner(u,t)*basis.subs( x == GiNaC::symbol("xi[0]"))
.subs( y == GiNaC::symbol("xi[1]")), d));
}
}
}
if ( order%2==1 )
{
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::lst tangent_vec = tangent(triangle, i);
GiNaC::ex Pk_space = bernstein(order-1, line, istr("a",i));
GiNaC::ex Pk = Pk_space.op(2);
GiNaC::ex Vt = inner(V, tangent_vec);
GiNaC::ex basis;
for (unsigned int j=0; j< Pk.nops(); j++)
{
basis = Pk.op(j);
GiNaC::ex integrand = Vt*basis;
GiNaC::ex dofi = line.integrate(integrand);
// dofs.insert(dofs.end(), dofi);
GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),2);
dofs.insert(dofs.end(), d);
GiNaC::ex eq = dofi == GiNaC::numeric(0);
equations.append(eq);
GiNaC::ex u = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]")));
GiNaC::ex t = GiNaC::matrix(2,1,GiNaC::lst(-GiNaC::symbol("normal_vec[1]"), GiNaC::symbol("normal_vec[0]")));
dof_repr.append(GiNaC::lst(inner(u,t)*basis.subs( x == GiNaC::symbol("xi[0]"))
.subs( y == GiNaC::symbol("xi[1]")), d));
}
}
}
// dofs related to the whole triangle
GiNaC::lst bernstein_polv;
if ( order > 0)
{
bernstein_polv = bernsteinv(2,order-1, triangle, "a");
GiNaC::ex basis_space = bernstein_polv.op(2);
for (unsigned int i=0; i< basis_space.nops(); i++)
{
GiNaC::lst basis = GiNaC::ex_to<GiNaC::lst>(basis_space.op(i));
GiNaC::ex integrand = inner(V, basis);
GiNaC::ex dofi = triangle.integrate(integrand);
GiNaC::ex eq = dofi == GiNaC::numeric(0);
equations.append(eq);
GiNaC::lst d = GiNaC::lst(GiNaC::lst(triangle.vertex(0), triangle.vertex(1), triangle.vertex(2)), i);
dofs.insert(dofs.end(), d);
}
}
variables = collapse(GiNaC::lst(V_vars.op(0),V_vars.op(1)));
// std::cout <<"no variables "<<variables.nops()<<std::endl;
// std::cout <<"no equations "<<equations.nops()<<std::endl;
GiNaC::matrix b; GiNaC::matrix A;
matrix_from_equations(equations, variables, A, b);
// std::cout <<" A "<<A.evalf()<<std::endl;
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(11+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 = V.op(0).subs(subs);
GiNaC::ex Nj2 = V.op(1).subs(subs);
Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2)));
b.let_op(11+i) = GiNaC::numeric(0);
}
}
}
}
| void SyFi::Robust::compute_basis_functions_old | ( | ) | [virtual] |
Definition at line 510 of file Robust.cpp.
References SyFi::bernstein(), SyFi::bernsteinv(), test_syfi::debug::c, SyFi::coeff(), SyFi::collapse(), SyFi::StandardFE::description, SyFi::div(), dof_repr, SyFi::StandardFE::dofs, SyFi::inner(), SyFi::Line::integrate(), SyFi::istr(), SyFi::Triangle::line(), SyFi::matrix_from_equations(), SyFi::normal(), SyFi::StandardFE::Ns, SyFi::StandardFE::order, SyFi::StandardFE::p, SyFi::pol2basisandcoeff(), SyFi::Line::repr(), SyFi::Polygon::str(), SyFi::t, SyFi::tangent(), SyFi::Polygon::vertex(), SyFi::x, and SyFi::y.
{
// remove previously computed basis functions and dofs
Ns.clear();
dofs.clear();
order = 3;
if ( p == NULL )
{
throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
}
if ( p->str().find("Line") != string::npos )
{
cout <<"Can not define the Robust element on a line"<<endl;
}
else if ( p->str().find("Triangle") != string::npos )
{
description = "Robust_2D";
Triangle& triangle = (Triangle&)(*p);
GiNaC::lst equations;
GiNaC::lst variables;
GiNaC::ex V_space = bernsteinv(2, order, triangle, "a");
GiNaC::ex V = V_space.op(0);
GiNaC::ex V_vars = V_space.op(1);
GiNaC::ex divV = div(V);
exmap basis2coeff = pol2basisandcoeff(divV);
exmap::iterator iter;
// div constraints:
for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
{
GiNaC::ex basis = (*iter).first;
GiNaC::ex coeff= (*iter).second;
if ( coeff != 0 && ( basis.degree(x) > 0 || basis.degree(y) > 0 ) )
{
equations.append( coeff == 0 );
}
}
// constraints on edges:
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::symbol s("s");
GiNaC::lst normal_vec = normal(triangle, i);
GiNaC::ex Vn = inner(V, normal_vec);
Vn = Vn.subs(line.repr(s).op(0)).subs(line.repr(s).op(1));
basis2coeff = pol2basisandcoeff(Vn,s);
for (iter = basis2coeff.begin(); iter != basis2coeff.end(); iter++)
{
GiNaC::ex basis = (*iter).first;
GiNaC::ex coeff= (*iter).second;
if ( coeff != 0 && basis.degree(s) > 1 )
{
equations.append( coeff == 0 );
}
}
}
// dofs related to the normal on the edges
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::lst normal_vec = normal(triangle, i);
GiNaC::ex Pk1_space = bernstein(1, line, istr("a",i));
GiNaC::ex Pk1 = Pk1_space.op(2);
GiNaC::ex Vn = inner(V, normal_vec);
GiNaC::ex basis;
for (unsigned int j=0; j< Pk1.nops(); j++)
{
basis = Pk1.op(j);
GiNaC::ex integrand = Vn*basis;
GiNaC::ex dofi = line.integrate(integrand);
// dofs.insert(dofs.end(), dofi);
GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),j);
dofs.insert(dofs.end(), d);
GiNaC::ex eq = dofi == GiNaC::numeric(0);
equations.append(eq);
GiNaC::ex u = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]")));
GiNaC::ex n = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("normal_vec[0]"), GiNaC::symbol("normal_vec[1]")));
dof_repr.append(GiNaC::lst(inner(u,n)*basis.subs( x == GiNaC::symbol("xi[0]"))
.subs( y == GiNaC::symbol("xi[1]")), d));
}
}
// dofs related to the tangent on the edges
for (int i=0; i< 3; i++)
{
Line line = triangle.line(i);
GiNaC::lst tangent_vec = tangent(triangle, i);
GiNaC::ex Pk_space = bernstein(0, line, istr("a",i));
GiNaC::ex Pk = Pk_space.op(2);
GiNaC::ex Vt = inner(V, tangent_vec);
GiNaC::ex basis;
for (unsigned int j=0; j< Pk.nops(); j++)
{
basis = Pk.op(j);
GiNaC::ex integrand = Vt*basis;
GiNaC::ex dofi = line.integrate(integrand);
// dofs.insert(dofs.end(), dofi);
GiNaC::lst d = GiNaC::lst(GiNaC::lst(line.vertex(0), line.vertex(1)),2);
dofs.insert(dofs.end(), d);
GiNaC::ex eq = dofi == GiNaC::numeric(0);
equations.append(eq);
GiNaC::ex u = GiNaC::matrix(2,1,GiNaC::lst(GiNaC::symbol("v[0]"), GiNaC::symbol("v[1]")));
GiNaC::ex t = GiNaC::matrix(2,1,GiNaC::lst(-GiNaC::symbol("normal_vec[1]"), GiNaC::symbol("normal_vec[0]")));
dof_repr.append(GiNaC::lst(inner(u,t)*basis.subs( x == GiNaC::symbol("xi[0]"))
.subs( y == GiNaC::symbol("xi[1]")), d));
}
}
variables = collapse(GiNaC::lst(V_vars.op(0),V_vars.op(1)));
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(11+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 = V.op(0).subs(subs);
GiNaC::ex Nj2 = V.op(1).subs(subs);
Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2)));
b.let_op(11+i) = GiNaC::numeric(0);
}
}
}
| GiNaC::lst SyFi::Robust::dof_repr |
Definition at line 30 of file Robust.h.
Referenced by compute_basis_functions(), and compute_basis_functions_old().
Definition at line 29 of file Robust.h.
Referenced by compute_basis_functions().
1.7.6.1