basis_symbolic.cc
Go to the documentation of this file.
#include "basis_symbolic.h"
#include "cln/cln.h"
using namespace std;
using namespace rheolef;
using namespace GiNaC;
basis_symbolic_nodal_on_geo::value_type
basis_symbolic_nodal_on_geo::eval (
const polynom_type& p,
const point_basic<ex>& xi,
size_type d) const
{
ex expr = p;
if (d > 0) expr = expr.subs(x == ex(xi[0]));
if (d > 1) expr = expr.subs(y == ex(xi[1]));
if (d > 2) expr = expr.subs(z == ex(xi[2]));
return expr;
}
matrix
basis_symbolic_nodal_on_geo::vandermonde_matrix (
const vector<ex>& p,
size_type d) const
{
unsigned int n = _node.size();
matrix vdm (n, n);
for (unsigned int i = 0; i < n; i++) {
const point_basic<ex>& xi = _node[i];
for (unsigned int j = 0; j < n; j++) {
vdm(i,j) = eval (p[j], xi, d);
}
}
return vdm;
}
/*
* Build the Lagrange basis, associated to nodes.
*
* input: nodes[n] and polynomial_basis[n]
* output: node_basis[n]
* such that node_basis[j](node[i]) = kronecker[i][j]
*
* algorithm: let:
* b_i = \sum_i a_{i,j} p_j
* where:
* p_j = polynomial basis [1, x, y, ..]
* b_i = basis associated to node :
* we want:
* b_i(x_k) = delta_{i,k}
* <=>
* a_{i,j} p_j(x_k) = \delta_{i,k}
* Let A = (a_{k,j})_{i,j} and c_{i,j} = (p_j(x_i))
* Then a_{i,j} c_{k,j} = delta_{i,k}
* <=>
* A = C^{-T}
*/
void
basis_symbolic_nodal_on_geo::make_node_basis()
{
assert_macro (_node.size() == _poly.size(),
"incompatible node set size (" << _node.size()
<< ") and polynomial basis size (" << _poly.size() << ").");
const size_type d = _hat_K.dimension();
const size_type n = size();
// Vandermonde matrix vdm(i,j) = pj(xi)
matrix vdm = vandermonde_matrix (_poly, d);
ex det_vdm = determinant(vdm);
if (det_vdm == 0) {
warning_macro("basis unisolvence failed on element `" << _hat_K.name() << "'");
for (size_type i = 0, n = _node.size(); i < n; i++) {
cerr << "node("<<i<<") = ";
_node[i].put (cerr, d);
cerr << endl;
}
for (size_type i = 0, n = _node.size(); i < n; i++) {
cerr << "poly("<<i<<") = " << _poly[i] << endl;
}
error_macro("basis unisolvence failed: unrecoverable.");
}
int det_vdm_exponent = floor(ex_to<numeric>(log(1.0*abs(det_vdm))).to_double()/log(10.0));
double det_vdm_mantisse = ex_to<numeric>(det_vdm/pow(ex(10.0),ex(det_vdm_exponent))).to_double();
warning_macro (_name<<"("<<_hat_K.name()<<"): determinant(vdm("<<n<<","<<n<<"))="<<ex_to<numeric>(det_vdm_mantisse).to_double() << "e" << det_vdm_exponent);
matrix inv_vdm = vdm.inverse();
// basis := trans(a)*poly
_basis.resize(n);
for (size_type i = 0; i < n; i++) {
polynom_type s = 0;
for (size_type j = 0; j < n; j++) {
s += inv_vdm(j,i) * _poly[j];
}
s = expand(s);
s = normal(s);
_basis[i] = s;
}
#ifdef VERY_VERBOSE
for (size_type i = 0; i < _basis.size(); i++) {
cerr << _hat_K.name() << ".basis("<<i<<") = " << _basis[i] << flush << endl;
}
#endif // VERY_VERBOSE
// check:
matrix vdm_l = vandermonde_matrix (_basis, d);
int ndigit10 = Digits;
numeric tol = ex_to<numeric>(pow(10.,-ndigit10/2.));
int status = 0;
for (size_type i = 0; i < n; i++) {
for (size_type j = 0; j < n; j++) {
if ((i == j && abs(vdm_l(i,j) - 1) > tol) ||
(i != j && abs(vdm_l(i,j)) > tol) ) {
error_macro ("Lagrange polynom check failed.");
}
}
}
// derivatives of the basis
_grad_basis.resize(n);
for (size_type i = 0; i < n; i++) {
if (d > 0) _grad_basis [i][0] = _basis[i].diff(x).expand().normal();
if (d > 1) _grad_basis [i][1] = _basis[i].diff(y).expand().normal();
if (d > 2) _grad_basis [i][2] = _basis[i].diff(z).expand().normal();
}
#ifdef VERY_VERBOSE
for (size_type i = 0; i < _basis.size(); i++) {
cerr << _hat_K.name() << ".grad_basis("<<i<<") = [" << flush;
for (size_type j = 0; j < d; j++) {
cerr << _grad_basis [i][j];
if (j != d-1) cerr << ", " << flush;
}
cerr << "]" << endl << flush;
}
#endif // VERY_VERBOSE
}
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
rheolef::point_basic
Definition: point.h:87
rheolef::normal
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation
Definition: field_expr_terminal.h:439
rheolef::basis_symbolic_nodal_on_geo::polynom_type
GiNaC::ex polynom_type
Definition: basis_symbolic.h:45
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
p
Definition: sphere.icc:25
rheolef::determinant
U determinant(const tensor_basic< U > &A, size_t d=3)
assert_macro
#define assert_macro(ok_condition, message)
Definition: dis_macros.h:113
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
basis_symbolic.h
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
rheolef::basis_symbolic_nodal_on_geo::size_type
std::vector< int >::size_type size_type
Definition: basis_symbolic.h:44
rheolef::std
Definition: vec_expr_v2.h:391
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
rheolef::basis_symbolic_nodal_on_geo::value_type
GiNaC::ex value_type
Definition: basis_symbolic.h:46