Rheolef  7.1
an efficient C++ finite element environment
basis_symbolic.cc
Go to the documentation of this file.
1 #include "basis_symbolic.h"
22 #include "cln/cln.h"
23 using namespace std;
24 using namespace rheolef;
25 using namespace GiNaC;
26 
28 basis_symbolic_nodal_on_geo::eval (
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 }
39 matrix
40 basis_symbolic_nodal_on_geo::vandermonde_matrix (
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 
79 void
80 basis_symbolic_nodal_on_geo::make_node_basis()
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 }
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
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
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
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
basis_symbolic.h
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:402
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