25 using namespace GiNaC;
28 basis_symbolic_nodal_on_geo::eval (
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]));
40 basis_symbolic_nodal_on_geo::vandermonde_matrix (
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);
80 basis_symbolic_nodal_on_geo::make_node_basis()
83 "incompatible node set size (" << _node.size()
84 <<
") and polynomial basis size (" << _poly.size() <<
").");
90 matrix vdm = vandermonde_matrix (_poly,
d);
93 warning_macro(
"basis unisolvence failed on element `" << _hat_K.name() <<
"'");
94 for (
size_type i = 0,
n = _node.size(); i <
n; i++) {
95 cerr <<
"node("<<i<<
") = ";
96 _node[i].put (cerr,
d);
99 for (
size_type i = 0,
n = _node.size(); i <
n; 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];
121 for (
size_type i = 0; i < _basis.size(); i++) {
122 cerr << _hat_K.name() <<
".basis("<<i<<
") = " << _basis[i] << flush << endl;
124 #endif // VERY_VERBOSE
126 matrix vdm_l = vandermonde_matrix (_basis,
d);
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) ) {
140 _grad_basis.resize(
n);
142 if (
d > 0) _grad_basis [i][0] = _basis[i].diff(x).expand().normal();
143 if (
d > 1) _grad_basis [i][1] = _basis[i].diff(y).expand().normal();
144 if (
d > 2) _grad_basis [i][2] = _basis[i].diff(z).expand().normal();
147 for (
size_type i = 0; i < _basis.size(); i++) {
148 cerr << _hat_K.name() <<
".grad_basis("<<i<<
") = [" << flush;
150 cerr << _grad_basis [i][j];
151 if (j !=
d-1) cerr <<
", " << flush;
153 cerr <<
"]" << endl << flush;
155 #endif // VERY_VERBOSE