34 if (
d > 0) expr = expr.subs(
x == ex(xi[0]));
35 if (
d > 1) expr = expr.subs(
y == ex(xi[1]));
36 if (
d > 2) expr = expr.subs(
z == ex(xi[2]));
44 unsigned int n =
_node.size();
47 for (
unsigned int i = 0; i < n; i++) {
50 for (
unsigned int j = 0; j < n; j++) {
51 vdm(i,j) =
eval (
p[j], xi,
d);
83 "incompatible node set size (" <<
_node.size()
84 <<
") and polynomial basis size (" <<
_poly.size() <<
").");
95 cerr <<
"node("<<i<<
") = ";
100 cerr <<
"poly("<<i<<
") = " <<
_poly[i] << endl;
102 error_macro(
"basis unisolvence failed: unrecoverable.");
104 int det_vdm_exponent = floor(ex_to<numeric>(log(1.0*abs(det_vdm))).to_double()/log(10.0));
105 double det_vdm_mantisse = ex_to<numeric>(det_vdm/
pow(ex(10.0),ex(det_vdm_exponent))).to_double();
106 warning_macro (
_name<<
"("<<
_hat_K.
name()<<
"): determinant(vdm("<<n<<
","<<n<<
"))="<<ex_to<numeric>(det_vdm_mantisse).to_double() <<
"e" << det_vdm_exponent);
107 matrix inv_vdm = vdm.inverse();
114 s += inv_vdm(j,i) *
_poly[j];
122 cerr <<
_hat_K.
name() <<
".basis("<<i<<
") = " <<
_basis[i] << flush << endl;
127 int ndigit10 = Digits;
129 numeric tol = ex_to<numeric>(
pow(10.,-ndigit10/2.));
133 if ((i == j && abs(vdm_l(i,j) - 1) > tol) ||
134 (i != j && abs(vdm_l(i,j)) > tol) ) {
148 cerr <<
_hat_K.
name() <<
".grad_basis("<<i<<
") = [" << flush;
151 if (j !=
d-1) cerr <<
", " << flush;
153 cerr <<
"]" << endl << flush;
std::vector< polynom_type > _poly
std::vector< point_basic< polynom_type > > _grad_basis
std::vector< point_basic< GiNaC::ex > > _node
std::vector< polynom_type > _basis
value_type eval(const polynom_type &p, const point_basic< polynom_type > &x, size_type d=3) const
std::vector< int >::size_type size_type
GiNaC::matrix vandermonde_matrix(const std::vector< polynom_type > &p, size_type d=3) const
size_type dimension() const
#define assert_macro(ok_condition, message)
#define error_macro(message)
#define warning_macro(message)
This file is part of Rheolef.
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
U determinant(const tensor_basic< U > &A, size_t d=3)
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation