Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
basis_symbolic.cc
Go to the documentation of this file.
1
21#include "basis_symbolic.h"
22#include "cln/cln.h"
23using namespace std;
24using namespace rheolef;
25using namespace GiNaC;
26
29 const polynom_type& p,
30 const point_basic<ex>& xi,
31 size_type d) const
32{
33 ex expr = p;
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]));
37 return expr;
38}
39matrix
41 const vector<ex>& p,
42 size_type d) const
43{
44 unsigned int n = _node.size();
45 matrix vdm (n, n);
46
47 for (unsigned int i = 0; i < n; i++) {
48
49 const point_basic<ex>& xi = _node[i];
50 for (unsigned int j = 0; j < n; j++) {
51 vdm(i,j) = eval (p[j], xi, d);
52 }
53 }
54 return vdm;
55}
56
57/*
58 * Build the Lagrange basis, associated to nodes.
59 *
60 * input: nodes[n] and polynomial_basis[n]
61 * output: node_basis[n]
62 * such that node_basis[j](node[i]) = kronecker[i][j]
63 *
64 * algorithm: let:
65 * b_i = \sum_i a_{i,j} p_j
66 * where:
67 * p_j = polynomial basis [1, x, y, ..]
68 * b_i = basis associated to node :
69 * we want:
70 * b_i(x_k) = delta_{i,k}
71 * <=>
72 * a_{i,j} p_j(x_k) = \delta_{i,k}
73 * Let A = (a_{k,j})_{i,j} and c_{i,j} = (p_j(x_i))
74 * Then a_{i,j} c_{k,j} = delta_{i,k}
75 * <=>
76 * A = C^{-T}
77 */
78
79void
81{
82 assert_macro (_node.size() == _poly.size(),
83 "incompatible node set size (" << _node.size()
84 << ") and polynomial basis size (" << _poly.size() << ").");
85
86 const size_type d = _hat_K.dimension();
87 const size_type n = size();
88
89 // Vandermonde matrix vdm(i,j) = pj(xi)
90 matrix vdm = vandermonde_matrix (_poly, d);
91 ex det_vdm = determinant(vdm);
92 if (det_vdm == 0) {
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);
97 cerr << endl;
98 }
99 for (size_type i = 0, n = _node.size(); i < n; i++) {
100 cerr << "poly("<<i<<") = " << _poly[i] << endl;
101 }
102 error_macro("basis unisolvence failed: unrecoverable.");
103 }
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();
108
109 // basis := trans(a)*poly
110 _basis.resize(n);
111 for (size_type i = 0; i < n; i++) {
112 polynom_type s = 0;
113 for (size_type j = 0; j < n; j++) {
114 s += inv_vdm(j,i) * _poly[j];
115 }
116 s = expand(s);
117 s = normal(s);
118 _basis[i] = s;
119 }
120#ifdef VERY_VERBOSE
121 for (size_type i = 0; i < _basis.size(); i++) {
122 cerr << _hat_K.name() << ".basis("<<i<<") = " << _basis[i] << flush << endl;
123 }
124#endif // VERY_VERBOSE
125 // check:
126 matrix vdm_l = vandermonde_matrix (_basis, d);
127 int ndigit10 = Digits;
128
129 numeric tol = ex_to<numeric>(pow(10.,-ndigit10/2.));
130 int status = 0;
131 for (size_type i = 0; i < n; i++) {
132 for (size_type j = 0; j < n; j++) {
133 if ((i == j && abs(vdm_l(i,j) - 1) > tol) ||
134 (i != j && abs(vdm_l(i,j)) > tol) ) {
135 error_macro ("Lagrange polynom check failed.");
136 }
137 }
138 }
139 // derivatives of the basis
140 _grad_basis.resize(n);
141 for (size_type i = 0; i < n; i++) {
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();
145 }
146#ifdef VERY_VERBOSE
147 for (size_type i = 0; i < _basis.size(); i++) {
148 cerr << _hat_K.name() << ".grad_basis("<<i<<") = [" << flush;
149 for (size_type j = 0; j < d; j++) {
150 cerr << _grad_basis [i][j];
151 if (j != d-1) cerr << ", " << flush;
152 }
153 cerr << "]" << endl << flush;
154 }
155#endif // VERY_VERBOSE
156}
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
#define assert_macro(ok_condition, message)
Definition dis_macros.h:113
#define error_macro(message)
Definition dis_macros.h:49
#define warning_macro(message)
Definition dis_macros.h:53
This file is part of Rheolef.
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
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
STL namespace.
Definition sphere.icc:25