|
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 () |
| def | __init__ |
| def | compute_basis_functions_old |
Public Attributes | |
| bool | pointwise |
| GiNaC::lst | dof_repr |
Static Public Attributes | |
| tuple | thisown = _swig_property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc='The membership flag') |
| tuple | pointwise = _swig_property(_SyFi.Robust_pointwise_get, _SyFi.Robust_pointwise_set) |
| tuple | dof_repr = _swig_property(_SyFi.Robust_dof_repr_get, _SyFi.Robust_dof_repr_set) |
Static Private Attributes | |
| __repr__ = _swig_repr | |
| __swig_destroy__ = _SyFi.delete_Robust | |
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] |
| def SyFi::Robust::__init__ | ( | self, | |
| args | |||
| ) |
__init__(SyFi::Robust self) -> Robust __init__(SyFi::Robust self, Polygon p, unsigned int order=0, bool pointwise=True) -> Robust __init__(SyFi::Robust self, Polygon p, unsigned int order=0) -> Robust __init__(SyFi::Robust self, Polygon p) -> Robust
Reimplemented from SyFi::StandardFE.
Definition at line 2495 of file SyFi.py.
02495 02496 def __init__(self, *args): 02497 """ 02498 __init__(SyFi::Robust self) -> Robust 02499 __init__(SyFi::Robust self, Polygon p, unsigned int order=0, bool pointwise=True) -> Robust 02500 __init__(SyFi::Robust self, Polygon p, unsigned int order=0) -> Robust 02501 __init__(SyFi::Robust self, Polygon p) -> Robust 02502 """ _SyFi.Robust_swiginit(self,_SyFi.new_Robust(*args))
| void SyFi::Robust::compute_basis_functions | ( | ) | [virtual] |
Reimplemented from SyFi::StandardFE.
Definition at line 41 of file Robust.cpp.
References SyFi.exmap::begin(), SyFi::bernstein(), SyFi::bernsteinv(), test_syfi::debug::c, SyFi::coeff(), SyFi::collapse(), SyFi::StandardFE::description, SyFi::div(), dof_repr, SyFi::StandardFE::dofs, SyFi.exmap::end(), SyFi::inner(), SyFi::Line::integrate(), SyFi::Triangle::integrate(), SyFi::interior_coordinates(), SyFi::istr(), SyFi.exmap::iterator(), 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.exmap::begin(), SyFi::bernstein(), SyFi::bernsteinv(), test_syfi::debug::c, SyFi::coeff(), SyFi::collapse(), SyFi::StandardFE::description, SyFi::div(), dof_repr, SyFi::StandardFE::dofs, SyFi.exmap::end(), SyFi::inner(), SyFi::Line::integrate(), SyFi::istr(), SyFi.exmap::iterator(), 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.
Referenced by compute_basis_functions_old().
{
// 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);
}
}
}
| def SyFi::Robust::compute_basis_functions_old | ( | self | ) |
compute_basis_functions_old(Robust self)
Definition at line 2504 of file SyFi.py.
References compute_basis_functions_old(), and SyFi.new_instancemethod.
02504 02505 def compute_basis_functions_old(self): 02506 """compute_basis_functions_old(Robust self)""" 02507 return _SyFi.Robust_compute_basis_functions_old(self) 02508 Robust.compute_basis_functions_old = new_instancemethod(_SyFi.Robust_compute_basis_functions_old,None,Robust)
SyFi::Robust::__repr__ = _swig_repr [static, private] |
Reimplemented from SyFi::StandardFE.
SyFi::Robust::__swig_destroy__ = _SyFi.delete_Robust [static, private] |
Reimplemented from SyFi::StandardFE.
| GiNaC::lst SyFi::Robust::dof_repr |
Definition at line 30 of file Robust.h.
Referenced by compute_basis_functions(), and compute_basis_functions_old().
tuple SyFi::Robust::dof_repr = _swig_property(_SyFi.Robust_dof_repr_get, _SyFi.Robust_dof_repr_set) [static] |
Definition at line 29 of file Robust.h.
Referenced by compute_basis_functions().
tuple SyFi::Robust::pointwise = _swig_property(_SyFi.Robust_pointwise_get, _SyFi.Robust_pointwise_set) [static] |
tuple SyFi::Robust::thisown = _swig_property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc='The membership flag') [static] |
Reimplemented from SyFi::StandardFE.
1.7.6.1