Rheolef  7.1
an efficient C++ finite element environment
basis_fem_RTk.cc
Go to the documentation of this file.
1 #include "basis_fem_RTk.h"
22 #include "equispaced.icc"
23 #include "warburton.icc"
24 #include "basis_ordering.icc"
25 #include "reference_element_aux.icc"
26 #include "eigen_util.h"
28 #include "piola_fem_hdiv.h"
29 
30 namespace rheolef {
31 using namespace std;
32 
33 // =========================================================================
34 // basis members
35 // =========================================================================
36 // allocators
37 // -------------------------------------------------------------------------
38 template<class T>
40 {
41 }
42 template<class T>
44  : basis_rep<T> (sopt),
45  _nnod_on_subgeo_remainder(),
46  _b_pre_kp1(),
47  _hat_node(),
48  _vdm(),
49  _inv_vdm(),
50  _quad(),
51  _bar_a(),
52  _bkm1_node_internal_d()
53 {
54 #ifdef TO_CLEAN
56 #endif // TO_CLEAN
58  // requieres a hierarchical pre basis => Dubiner or monomial
59  // - monomial is ill conditionned, but easier to decompose
60  // - dubiner: TODO project the tilde_psi basis on it
61  // note: the internal dof basis is independent
62  _b_pre_kp1 = basis_raw_basic<T> ("M"+itos(k+1));
64 
65  // piola FEM transformation:
66  typedef piola_fem_hdiv<T> piola_fem_type;
67  base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_type));
68 }
69 template<class T>
70 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
72 {
73  base::_initialize_data_guard (hat_K);
74  return _hat_node [hat_K.variant()];
75 }
76 template<class T>
77 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
79 {
80  base::_initialize_data_guard (hat_K);
81  return _vdm [hat_K.variant()];
82 }
83 template<class T>
84 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
86 {
87  base::_initialize_data_guard (hat_K);
88  return _inv_vdm [hat_K.variant()];
89 }
90 template<class T>
91 void
93 {
94  size_type k = degree()-1;
95  for (size_type map_d = 0; map_d < 4; ++map_d) {
96  base::_ndof_on_subgeo [map_d].fill (0);
97  base::_nnod_on_subgeo [map_d].fill (0);
98  }
99  // --------------------------------------------------------------------------
100  // 1. compute ndof and dof pointers:
101  // --------------------------------------------------------------------------
102  // see BreFor-1991 p. 116 for K=t,T & p. 119 for K=q,H ; BerDur-2013 for K=P
103  base::_ndof_on_subgeo [0][reference_element::p] = 1;
104  base::_ndof_on_subgeo [1][reference_element::p] = 1;
105  base::_ndof_on_subgeo [1][reference_element::e] = k; // RT0 is P1 on 1d with K=e
106  base::_ndof_on_subgeo [2][reference_element::e] = k+1;
107  base::_ndof_on_subgeo [2][reference_element::t] = k*(k+1);
108  base::_ndof_on_subgeo [2][reference_element::q] = 2*k*(k+1); // Rav-1981, p. 23
109  base::_ndof_on_subgeo [3][reference_element::t] = (k+1)*(k+2)/2;
110  base::_ndof_on_subgeo [3][reference_element::q] = sqr(k+1);
111  base::_ndof_on_subgeo [3][reference_element::T] = 3*k*(k+1)*(k+2)/6;
112  base::_ndof_on_subgeo [3][reference_element::P] = sqr(k)*(k+1)/2; // see BerDur-2013
113  base::_ndof_on_subgeo [3][reference_element::H] = 3*k*sqr(k+1);
114  base::_helper_discontinuous_ndof_on_subgeo_inplace_change (base::is_continuous(), base::_ndof_on_subgeo);
115  base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_ndof_on_subgeo, base::_first_idof_by_dimension);
116  // --------------------------------------------------------------------------
117  // 2. nodes for dofs : interpolate v.n on sides + internal quadrature nodes
118  // --------------------------------------------------------------------------
119  // nodes coincides with dofs on sides : Lagrange nodes
120  // interior nodes are quadrature node
121  quadrature_option qopt;
122  qopt.set_family (quadrature_option::gauss);
123  // order = 2*(k-1)) for P_{k-1} internal modes for simplicial elts hat_K=etT
124  // = 2*k for P_{k-1,k,k} internal modes for non-simplicial elts hat_K=qHP
125  qopt.set_order (2*(k+2)); // takes the max to have an uniform quadrature for any cases TODO: too much ? min-adjust=2*k to converge
126  _quad = quadrature<T> (qopt);
127  base::_nnod_on_subgeo [0][reference_element::p] = 1;
128  base::_nnod_on_subgeo [1][reference_element::p] = 1;
129  base::_nnod_on_subgeo [1][reference_element::e] = k == 0 ? 0 : _quad.size (reference_element::e);
130  base::_nnod_on_subgeo [2][reference_element::e] = k+1;
131  base::_nnod_on_subgeo [2][reference_element::t] = k*(k+1)/2; // k == 0 ? 0 : _quad.size (reference_element::t);
132  base::_nnod_on_subgeo [2][reference_element::q] = k == 0 ? 0 : _quad.size (reference_element::q);
133  base::_nnod_on_subgeo [3][reference_element::t] = (k+1)*(k+2)/2;
134  base::_nnod_on_subgeo [3][reference_element::q] = sqr(k+1);
135  base::_nnod_on_subgeo [3][reference_element::T] = k == 0 ? 0 : _quad.size (reference_element::T);
136  base::_nnod_on_subgeo [3][reference_element::P] = k == 0 ? 0 : _quad.size (reference_element::P);
137  base::_nnod_on_subgeo [3][reference_element::H] = k == 0 ? 0 : _quad.size (reference_element::H);
138  _nnod_on_subgeo_remainder = base::_nnod_on_subgeo; // used for node table indexing
139  base::_helper_discontinuous_ndof_on_subgeo_inplace_change (base::is_continuous(), base::_nnod_on_subgeo);
140  base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_nnod_on_subgeo, base::_first_inod_by_dimension);
141 }
142 template<class T>
143 void
145 {
146  if (hat_K.dimension() == 0) return; // nothing to do
147 
148  // abbrevs
149  size_type d = hat_K.dimension();
150  size_type k = degree()-1;
151  size_type variant = hat_K.variant();
152 
153  // --------------------------------------------------------------------------
154  // 1. nodes for dofs : interpolate v.n on sides + internal quadrature nodes
155  // --------------------------------------------------------------------------
156  // 1.1. insert nodes
157  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_node = _hat_node [hat_K.variant()];
158  hat_node.resize (base::_first_inod_by_dimension [variant][d+1]);
159  size_type loc_nnod = hat_node.size();
160  // 1.2. nodes on sides
161  size_type loc_inod = 0;
162  for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(d-1); loc_isid < loc_nsid; ++loc_isid) {
163  reference_element hat_S = hat_K.subgeo (d-1, loc_isid);
164  size_type sid_variant = hat_S.variant();
165  size_type loc_nsidvert = hat_S.n_vertex();
166  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> hat_node_sid;
167  const size_type internal = 1;
168  switch (base::_sopt.get_node()) {
170 #ifdef TODO
171  // TODO: pointset warburton could support the "internal=1" option
172  pointset_lagrange_warburton (sid_variant, k, hat_node_sid, internal); break;
173 #endif // TODO
175  pointset_lagrange_equispaced (sid_variant, k, hat_node_sid, internal); break;
176  default:
177  error_macro ("unsupported node set: " << base::_sopt.get_node_name());
178  }
179  point loc_vertex [4]; // restricted to point_basic<Float> in reference_element::vertex()
180  for (size_type loc_jsidvert = 0; loc_jsidvert < loc_nsidvert; ++loc_jsidvert) {
181  size_type loc_jvertex = hat_K.subgeo_local_vertex (d-1, loc_isid, loc_jsidvert);
182  loc_vertex[loc_jsidvert] = hat_K.vertex (loc_jvertex);
183  }
184  for (size_type loc_inod_sid = 0, loc_nnod_sid = hat_node_sid.size(); loc_inod_sid < loc_nnod_sid; ++loc_inod_sid) {
185  check_macro (loc_inod < size_t(hat_node.size()), "invalid loc_inod");
186  T xi0 = hat_node_sid [loc_inod_sid][0],
187  xi1 = hat_node_sid [loc_inod_sid][1];
188  if (loc_nsidvert == 4) { // map from [-1,1]^2 to [0,1]^2
189  xi0 = (1+xi0)/2;
190  xi1 = (1+xi1)/2;
191  }
192  for (size_type alpha = 0; alpha < d; ++alpha) {
193  hat_node[loc_inod][alpha] = loc_vertex [0][alpha];
194  if (d >= 2) {
195  hat_node[loc_inod][alpha] += xi0*(loc_vertex[1][alpha] - loc_vertex[0][alpha]);
196  }
197  if (d == 3) {
198  hat_node[loc_inod][alpha] += xi1*(loc_vertex[loc_nsidvert-1][alpha] - loc_vertex[0][alpha]);
199  }
200  }
201  ++loc_inod;
202  }
203  }
204  // 1.3. insert internal quadrature nodes for integrating exactly P(k-1) in hat_K, when k>0
205  if (k > 0) {
206  if (variant == reference_element::t) {
207  // experimental for triangle: internal interpolation without quadrature
208  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> hat_node_internal;
209  const size_type internal = 1;
210  switch (base::_sopt.get_node()) {
212 #ifdef TODO
213  // TODO: pointset warburton could support the "internal=1" option
214  pointset_lagrange_warburton (variant, k-1, hat_node_internal, internal); break;
215 #endif // TODO
217  pointset_lagrange_equispaced (variant, k-1, hat_node_internal, internal); break;
218  default:
219  error_macro ("unsupported node set: " << base::_sopt.get_node_name());
220  }
221  trace_macro ("hat_node_internal.size="<<hat_node_internal.size());
222  for (size_type loc_inod_int = 0; loc_inod_int < size_type(hat_node_internal.size()); ++loc_inod_int) {
223  trace_macro ("hat_node_internal["<<loc_inod_int<<"]="<<hat_node_internal[loc_inod_int]);
224  hat_node[loc_inod] = hat_node_internal [loc_inod_int];
225  ++loc_inod;
226  }
227  } else { // variant == TPH : still quadrature
228  for (typename quadrature<T>::const_iterator iter_q = _quad.begin(hat_K),
229  last_q = _quad.end(hat_K); iter_q != last_q; iter_q++) {
230  hat_node[loc_inod] = (*iter_q).x;
231  ++loc_inod;
232  }
233  }
234  }
235  check_macro (loc_inod == loc_nnod, "invalid node count: loc_inod="<<loc_inod<<" and loc_nnod="<<loc_nnod);
236 
237  // --------------------------------------------------------------------------
238  // 2. build a transformation tilde_A for evaluating polynomials from (b_pre_{k+1})^d
239  // from a pre-basis:
240  // tilde_psi = {(p,0), (0,p), p in Bkm1} cup {(x0*p,x1*p), p in bar_Bk}
241  // --------------------------------------------------------------------------
242  // 2.1. compute ndof per side, for each side (prism have different sides)
243  size_type dim_Pkm1 = (k == 0) ? 0 : reference_element::n_node(hat_K.variant(), k-1);
244  size_type dim_Pk = reference_element::n_node (hat_K.variant(), k);
245  size_type dim_Pkp1 = _b_pre_kp1.ndof (hat_K);
246  size_type loc_ndof_sid_tot = 0;
247  for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(d-1); loc_isid < loc_nsid; ++loc_isid) {
248  reference_element hat_S = hat_K.subgeo (d-1, loc_isid);
249  // dofs on sides are Lagrange dofs, associated to nodes:
250  loc_ndof_sid_tot += _nnod_on_subgeo_remainder [d][hat_S.variant()];
251  }
252  size_type loc_ndof = base::_first_idof_by_dimension [variant][d+1];
253  size_type loc_nnod_sid_tot = loc_ndof_sid_tot; // dofs on sides are Lagrange one, associated to nodes
254 
255 #ifdef TO_CLEAN
256  // check size
257  warning_macro ("first_inod(hat_K,d) ="<<base::first_inod(hat_K,d));
258  warning_macro ("first_inod(hat_K,d+1)="<<base::first_inod(hat_K,d+1));
259  warning_macro ("dim(P"<<int(k-1)<<")="<<dim_Pkm1);
260  warning_macro ("dim(P"<<k<<")="<<dim_Pk);
261  warning_macro ("dim(P"<<k+1<<")="<<dim_Pkp1);
262  warning_macro ("loc_ndof_sid_tot="<<loc_ndof_sid_tot);
263  warning_macro ("size(hat_K)="<<base::ndof (hat_K));
264  warning_macro ("first_idof(hat_K,d-1)="<<base::first_idof(hat_K,d-1));
265  warning_macro ("first_idof(hat_K,d) ="<<base::first_idof(hat_K,d));
266  warning_macro ("first_idof(hat_K,d+1)="<<base::first_idof(hat_K,d+1));
267  warning_macro ("dim(RT"<<k<<")="<<loc_ndof);
268  warning_macro ("dim((P"<<k+1<<")^"<<d<<")="<<d*dim_Pkp1);
269  warning_macro ("a_tilde(loc_ndof)="<<loc_ndof<<",d*dim(Pkp1)="<<d*dim_Pkp1<<")");
270 #endif // TO_CLEAN
271 
272  check_macro (loc_ndof == base::ndof (hat_K), "invalid ndof="<<loc_ndof<<" != "<< base::ndof (hat_K));
273  //
274  // 2.2. decompose the tilde_psi RTk basis on the basis of (P_{k+1})^d
275  // note: explicit easy decomposition when using monomial basis for (P_{k+1})^d
276  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> tilde_a (loc_ndof, d*dim_Pkp1);
277  tilde_a.fill (0);
278  // 2.2.1 homogeneous polynoms:
279  // (q,0) and (0,r), q,r in basis(Pk)
280  // leads to a diag matrix block
281  // since basis(Pk) subset basis(P_{k+1}) and the basis is hierarchical
282  trace_macro ("basis(1): (q,0) and (0,r), q,r in basis(Pk)...");
283  for (size_type loc_comp_idof = 0; loc_comp_idof < dim_Pk; ++loc_comp_idof) {
284  for (size_type alpha = 0; alpha < d; ++alpha) {
285  size_type loc_idof = d*loc_comp_idof + alpha;
286  size_type loc_jpkp1 = d*loc_comp_idof + alpha;
287  trace_macro ("loc_idof="<<loc_idof<<", loc_jpkp1="<<loc_jpkp1);
288  tilde_a (loc_idof, loc_jpkp1) = 1;
289  }
290  }
291  // 2.2.2 non-homogeneous polynoms
292  trace_macro ("basis(2): (x0*p,x1*p), p in basis(P_k)\basis(P_{k-1}) [exactly k degree]");
293  std::vector<size_type> inod2ideg_kp1;
294  build_inod2ideg (hat_K, k+1, inod2ideg_kp1);
295  switch (hat_K.variant()) {
296  case reference_element::e:
297  case reference_element::t:
298  case reference_element::T: {
299  // triangle : (x0*p,x1*p), p in basis(P_k)\basis(P_{k-1}) [exactly k degree]
300  // =(b_iO,b_i1) in P_{k+1} when b is the hierarchical monomial basis
301  // and i0=ilat2ideg([i+1,k-i ]), i=0..k
302  // i1=ilat2ideg([i ,k-i+1])
303  std::vector<point_basic<size_type> > power_index;
304  make_power_indexes_sorted_by_degrees (hat_K, k, power_index);
305  for (size_type loc_ideg = dim_Pkm1, loc_ndeg = power_index.size(); loc_ideg < loc_ndeg; ++loc_ideg) {
306  for (size_type alpha = 0; alpha < d; ++alpha) {
307  point_basic<size_type> ilat = power_index [loc_ideg];
308  size_type loc_idof = d*dim_Pk + (loc_ideg - dim_Pkm1);
309  ilat [alpha]++; // x_alpha*monomial(i,k-i) = x^{i+1}*y^{k-i} when alpha=0, etc
310  size_type loc_inod_kp1 = ilat2loc_inod (hat_K, k+1, ilat);
311  size_type loc_ideg_kp1 = inod2ideg_kp1 [loc_inod_kp1];
312  size_type loc_jpkp1d = d*loc_ideg_kp1 + alpha;
313  trace_macro ("loc_idof="<<loc_idof<<", loc_jpkp1d="<<loc_jpkp1d);
314  tilde_a (loc_idof, loc_jpkp1d) = 1;
315  }
316  }
317  break;
318  }
320  case reference_element::H: {
321  // quadrangle : (x0*p, 0)
322  // ( ,x1*q), p,q in basis(P_k)\basis(P_{k-1}) [exactly k degree]
323  // =(b_i0, 0)
324  // ( 0,b_i1), b_i0, b_i1 in P_{k+1} when b is the hierarchical monomial basis
325  // and i0=ilat2ideg([i+1,j ]), i=0..k
326  // i1=ilat2ideg([i ,j+1])
327  size_type sub_variant = (d == 2) ? reference_element::e : reference_element::q;
328  reference_element hat_sub (sub_variant); // q=e*e ; H = e*q
329  std::vector<point_basic<size_type> > power_index_sub;
330  make_power_indexes_sorted_by_degrees (hat_sub, k, power_index_sub);
331  for (size_type loc_ideg = 0, loc_ndeg = power_index_sub.size(); loc_ideg < loc_ndeg; ++loc_ideg) {
332  for (size_type alpha = 0; alpha < d; ++alpha) {
333  point_basic<size_type> ilat_sub = power_index_sub [loc_ideg];
334  size_type loc_idof = d*dim_Pk + d*loc_ideg + alpha;
335  point_basic<size_type> ilat (0,0,0);
336  ilat [alpha] = k+1; // (x^{k+1}*p(y), 0) & (0, p(x)*y^{k+1}) with p in Pk(sub)
337  ilat [(alpha+1)%d] = ilat_sub[0];
338  if (d == 3) ilat [(alpha+2)%d] = ilat_sub[1];
339  size_type loc_inod_kp1 = ilat2loc_inod (hat_K, k+1, ilat);
340  size_type loc_ideg_kp1 = inod2ideg_kp1 [loc_inod_kp1];
341  size_type loc_jpkp1d = d*loc_ideg_kp1 + alpha;
342  trace_macro ("loc_idof="<<loc_idof<<", loc_jpkp1d="<<loc_jpkp1d);
343  tilde_a (loc_idof, loc_jpkp1d) = 1;
344  }
345  }
346  break;
347  }
349  default: error_macro ("unexpected element `"<<hat_K.name()<<"' (HINT: see BerDur-2013)");
350  }
351 #undef DEBUG_RTK // print matrix
352 #ifdef DEBUG_RTK
353  cout << setprecision(std::numeric_limits<T>::digits10)
354  << "tilde_a=[" << endl << tilde_a <<"]"<<endl
355  << "[u,s,vt]=svd(tilde_a)" << endl
356  << "id=eye(size(s))" << endl
357  << "a=u*id*vt" << endl
358  ;
359 #endif // DEBUG_RTK
360  // --------------------------------------------------------------------------
361  // 3. build a transformation A for evaluating polynomials from (b_pre_{k+1})^d
362  // for a raw-basis:
363  // psi = A*b_pre
364  // --------------------------------------------------------------------------
365  // get psi raw RTk basis by othogonalizing the tilde_psi RTk basis
366  // note: use SVD
367  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
368  u (tilde_a.rows(), tilde_a.rows()),
369  vt (tilde_a.cols(), tilde_a.cols());
370  Eigen::Matrix<T,Eigen::Dynamic,1> s (tilde_a.rows());
371  Eigen::JacobiSVD<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >
372  svd (tilde_a, Eigen::ComputeFullU | Eigen::ComputeFullV);
373  // SVD: tilde_a = u*s1*trans(vt) TODO check with eigen u/ut v/vt !
374  s = svd.singularValues();
375  u = svd.matrixU();
376  vt = svd.matrixV();
377  size_type rank_s = 0;
379  for (size_type loc_idof = 0; loc_idof < size_type(s.size()); ++loc_idof) {
380  rank_s += (abs(s[loc_idof]) > eps) ? 1 : 0;
381  }
382  check_macro (rank_s == loc_ndof,
383  "invalid polynomial space dimension = " << rank_s << " < loc_ndof = " << loc_ndof);
384  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> id (tilde_a.rows(), tilde_a.cols());
385  id.setIdentity();
386  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> a = id*vt.transpose();
387 
388 #ifdef DEBUG_RTK
389  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> s1 (tilde_a.rows(), tilde_a.cols());
390  s1.fill(0);
391  for (size_type iloc = 0; iloc < size_t(tilde_a.rows()); iloc++) {
392  s1(iloc,iloc) = s(iloc);
393  }
394  cout << "s1 = ["<< endl << s1 <<"];" << endl
395  << "u1 = ["<< endl << u <<"];"<<endl
396  << "vt1 = ["<< endl << vt <<"];"<<endl
397  << "id1 = ["<< endl << id <<"];"<<endl
398  << "a1 = ["<< endl << a <<"]"<<endl
399  << "err=norm(tilde_a-u*s*vt')"<<endl
400  << "err1=norm(tilde_a-u1*s1*vt1')"<<endl
401  << "err_a=abs(a-a1)"<<endl
402  << "err_svd=max(max(abs(a-a1)))"<<endl
403  << "err_u=max(max(abs(u-u1)))"<<endl
404  << "err_v=max(max(abs(vt-vt1)))"<<endl
405  << "err_s=max(diag(s1)-diag(s))"<<endl
406  ;
407  T err_svd = (tilde_a - u*s1*vt.transpose()).norm();
408  cout << "err_svd = " << err_svd << endl;
409 #endif // DEBUG_RTK
410 
411  // -------------------------------------------------------------------------------
412  // 4. build a transformation bar_A for evaluating polynomials from (b_pre_{k+1})^d
413  // for the dof-basis:
414  // phi = A*b_pre
415  // -------------------------------------------------------------------------------
416  // 4.1. defines the quadrature nodes, for integral dofs
417  // and evaluate the basis of P_{k-1} on it
418  std::array<
419  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
420  ,3>&
421  bkm1_node_internal_d = _bkm1_node_internal_d [variant];
422  if (k > 0) {
423  size_type loc_nnod_int = base::nnod_on_subgeo (d,hat_K.variant());
424  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> hat_node_internal (loc_nnod_int);
425  for (size_type loc_inod_int = 0; loc_inod_int < loc_nnod_int; ++loc_inod_int) {
426  size_type loc_inod = loc_nnod_sid_tot + loc_inod_int;
427  hat_node_internal [loc_inod_int] = hat_node [loc_inod];
428  }
429  string basis_dof_name;
430  switch (base::_sopt.get_raw_polynomial()) {
431  case basis_option::monomial: basis_dof_name = "M"; break;
432  case basis_option::dubiner: basis_dof_name = "D"; break;
433  case basis_option::bernstein: basis_dof_name = "B"; break;
434  default: error_macro ("unsupported raw polynomial basis `"<<base::_sopt.get_raw_polynomial_name()<<"'");
435  }
436 trace_macro("basis_dof_name="<<basis_dof_name);
437  // bkm1_node_internal(i,j) = bj(xi)
438  switch (variant) {
439  case reference_element::e:
440  case reference_element::t:
441  case reference_element::T: {
442  // int_K b*v dx with v=[p,0],[0,p] and p in P_{k-1}
443  basis_raw_basic<T> b_dof_km1 (basis_dof_name+itos(k-1));
444  details::basis_on_pointset_evaluate (b_dof_km1, hat_K, hat_node_internal, bkm1_node_internal_d[0]);
445  for (size_type alpha = 1; alpha < d; ++alpha) {
446  bkm1_node_internal_d [alpha] = bkm1_node_internal_d [0];
447  }
448  break;
449  }
451  case reference_element::H: {
452  // int_K b*v dx with v=[p,0],[0,q] and p in P_{k-1,k}, q in P_{k,k-1}
453  basis_raw_basic<T> b_dof_km1 (basis_dof_name+itos(k-1)),
454  b_dof_k (basis_dof_name+itos(k));
455  std::array<Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>,3> hat_node_internal_comp;
456  std::array<
457  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
458  ,3>
459  bkm1_node_internal_comp,
460  bk_node_internal_comp;
462  for (size_type alpha = 0; alpha < d; ++alpha) {
463  hat_node_internal_comp [alpha].resize (hat_node_internal.size());
464  for (size_type loc_inod_int = 0, loc_nnod_int = hat_node_internal.size(); loc_inod_int < loc_nnod_int; ++loc_inod_int) {
465  // resize from q=[-1,1]^2 to e=[0,1]
466  hat_node_internal_comp [alpha][loc_inod_int][0] = (1+hat_node_internal [loc_inod_int][alpha])/2;
467  }
468  details::basis_on_pointset_evaluate (b_dof_km1, hat_e, hat_node_internal_comp[alpha], bkm1_node_internal_comp[alpha]);
469  details::basis_on_pointset_evaluate (b_dof_k, hat_e, hat_node_internal_comp[alpha], bk_node_internal_comp[alpha]);
470  }
471  size_type loc_ndof_int = bkm1_node_internal_comp[0].cols()*pow(bk_node_internal_comp[0].cols(),d-1);
472  for (size_type alpha = 0; alpha < d; ++alpha) {
473  size_type alpha2 = (alpha+1)%d;
474  size_type alpha3 = (alpha+2)%d;
475  bkm1_node_internal_d [alpha].resize (hat_node_internal.size(), loc_ndof_int);
476  for (size_type loc_inod_int = 0, loc_nnod_int = hat_node_internal.size(); loc_inod_int < loc_nnod_int; ++loc_inod_int) {
477  size_type loc_idof_int = 0;
478  if (variant == reference_element::q) {
479  for (size_type i = 0; i < size_type(bkm1_node_internal_comp [alpha].cols()); ++i) {
480  for (size_type j = 0; j < size_type( bk_node_internal_comp [alpha].cols()); ++j) {
481  bkm1_node_internal_d [alpha ] (loc_inod_int,loc_idof_int)
482  = bkm1_node_internal_comp [alpha ] (loc_inod_int,i)
483  * bk_node_internal_comp [alpha2] (loc_inod_int,j);
484  ++loc_idof_int;
485  }}
486  } else {
487  for (size_type i = 0; i < size_type(bkm1_node_internal_comp [alpha].cols()); ++i) {
488  for (size_type j = 0; j < size_type( bk_node_internal_comp [alpha].cols()); ++j) {
489  for (size_type k = 0; k < size_type( bk_node_internal_comp [alpha].cols()); ++k) {
490  bkm1_node_internal_d [alpha ] (loc_inod_int,loc_idof_int)
491  = bkm1_node_internal_comp [alpha ] (loc_inod_int,i)
492  * bk_node_internal_comp [alpha2] (loc_inod_int,j)
493  * bk_node_internal_comp [alpha3] (loc_inod_int,k);
494  ++loc_idof_int;
495  }}}
496  }
497  }
498  }
499  break;
500  }
502  default: error_macro ("unexpected element `"<<hat_K.name()<<"'");
503  }
504  }
505  // -----------------------------------------
506  // 4.2. compute basis of (P_{k+1})^d at all nodes
507  // -----------------------------------------
508  size_type loc_n_bkp1 = _b_pre_kp1.ndof (hat_K);
509  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> bkp1_node (loc_nnod, loc_n_bkp1);
510  details::basis_on_pointset_evaluate (_b_pre_kp1, hat_K, hat_node, bkp1_node); // bkp1_node(i,j) = bj(xi)
511  // vector expansion: bkp1d_node(i,d*j+alpha) = [bj(xi),0,0], [0,bj(xi),0], etc
512  // -> compute the dofs for all this expansion
513  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> bkp1d_dof (loc_ndof, d*loc_n_bkp1); // dof_i(bjd)
514  bkp1d_dof.fill (std::numeric_limits<T>::max());
515  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> bkp1d_j_node (loc_nnod);
516  Eigen::Matrix<T,Eigen::Dynamic,1> bkp1d_j_dof (loc_ndof);
517  for (size_type loc_j_bkp1 = 0; loc_j_bkp1 < loc_n_bkp1; ++loc_j_bkp1) {
518  for (size_type alpha = 0; alpha < d; ++alpha) {
519  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
520  bkp1d_j_node [loc_inod] = point_basic<T>(0,0,0);
521  bkp1d_j_node [loc_inod][alpha] = bkp1_node(loc_inod,loc_j_bkp1);
522  }
523  bkp1d_j_dof.fill (std::numeric_limits<T>::max());
524  _compute_dofs (hat_K, bkp1d_j_node, bkp1d_j_dof);
525  size_type loc_j_bkp1d = d*loc_j_bkp1 + alpha;
526  check_macro (bkp1d_dof.rows() == bkp1d_j_dof.size(), "invalid sizes");
527  bkp1d_dof.col (loc_j_bkp1d) = bkp1d_j_dof;
528  }
529  }
530  // -----------------------------------------
531  // 4.3. vdm
532  // -----------------------------------------
533  // VDM(i,j) = dof_i(phi_j)
534  // = phi_j.n(xi) for side dofs
535  // int phi_j(x)*p(x)*dx, for all p in P(k-1)
536  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& vdm = _vdm [hat_K.variant()];
537  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& inv_vdm = _inv_vdm[hat_K.variant()];
538  vdm = bkp1d_dof*a.transpose(); // = trans(a*trans(bkp1d_dof));
539  bool invert_ok = invert(vdm, inv_vdm);
540 #ifdef TO_CLEAN
541  if (!invert_ok) {
542  warning_macro ("unisolvence failed for " << base::name() <<"(" << hat_K.name() << ") basis");
543  }
544 #endif // TO_CLEAN
545  check_macro (invert_ok,
546  "unisolvence failed for " << base::name() <<"(" << hat_K.name() << ") basis");
547  // -----------------------------------------
548  // 4.4. final composition matrix: bar_a = trans(inv_vdm)*a
549  // -----------------------------------------
550  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& bar_a = _bar_a [hat_K.variant()];
551  bar_a = inv_vdm.transpose()*a;
552 #ifdef DEBUG_RTK
553  cout << "bkp1d_dof=[" << endl << bkp1d_dof <<"]"<<endl
554  << "vdm=[" << endl << vdm <<"]"<<endl
555  << "det_vdm=" << vdm.determinant() <<endl
556  << "cond_vdm=" << cond(vdm) <<endl
557  << "bar_a1=inv(vdm)'*a;"<<endl
558  << "bar_a=[" << endl << bar_a <<"]"<<endl;
559 #endif // DEBUG_RTK
560 }
561 // ----------------------------------------------------------------------------
562 // evaluation of all basis functions at hat_x:
563 // ----------------------------------------------------------------------------
564 template<class T>
565 void
567  reference_element hat_K,
568  const point_basic<T>& hat_x,
569  Eigen::Matrix<value_type,Eigen::Dynamic,1>& value) const
570 {
571  base::_initialize_data_guard (hat_K);
572  size_type d = hat_K.dimension();
573  size_type loc_ndof = base::ndof (hat_K);
574  //
575  // 1) evaluate the basis of P_{k+1}(hat_K) at hat_x : bkp1(x)
576  //
577  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& bar_a = _bar_a [hat_K.variant()];
578  Eigen::Matrix<T,Eigen::Dynamic,1> bkp1;
579  _b_pre_kp1.evaluate (hat_K, hat_x, bkp1);
580  value.resize (loc_ndof);
581  //
582  // 2) evaluate the basis of RTk at hat_x: phi(x)
583  // [phi(x)] = [A2]*[bkp1d(x)]
584  // where bkp1d = basis of (P_{k+1}(hat_K))^d by vectorization
585  // nb
586  // ---
587  // \.
588  // phi_{i,a}(hat_x) = / A(i,[j,a])*b_j(hat_x)
589  // ---
590  // j=0
591  for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
592  value[loc_idof] = point_basic<T>(0,0,0);
593  for (size_type loc_jdof_bkp1 = 0, loc_ndof_bkp1 = bkp1.size(); loc_jdof_bkp1 < loc_ndof_bkp1; ++loc_jdof_bkp1) {
594  for (size_type alpha = 0; alpha < d; ++alpha) {
595  size_type loc_jdof_bkp1d = d*loc_jdof_bkp1 + alpha;
596  value[loc_idof][alpha] += bar_a(loc_idof,loc_jdof_bkp1d)*bkp1[loc_jdof_bkp1];
597  }
598  }
599  }
600 }
601 template<class T>
602 void
604  reference_element hat_K,
605  const point_basic<T>& hat_x,
606  Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& value) const
607 {
608  base::_initialize_data_guard (hat_K);
609  size_type d = hat_K.dimension();
610  size_type loc_ndof = base::ndof (hat_K);
611  //
612  // 1) evaluate the grad basis of P_{k+1}(hat_K) at hat_x : grad bkp1(x)
613  //
614  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& bar_a = _bar_a [hat_K.variant()];
615  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> grad_bkp1;
616  _b_pre_kp1.grad_evaluate (hat_K, hat_x, grad_bkp1);
617  value.resize (loc_ndof);
618  //
619  // 2) evaluate the grad basis of RTk at hat_x: grad phi(x)
620  // [grad phi(x)] = [A2]*[grad bkp1d(x)]
621  // where grad bkp1d = grad basis of (P_{k+1}(hat_K))^d by vectorization
622  //
623  // nb
624  // ---
625  // d phi_{i,a} \ d b_j
626  // ----------(hat_x) = / A(i,[j,a])*--------(hat_x)
627  // d hat_x_b --- d hat_x_b
628  // j=0
629  for (size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
630  value[loc_idof] = tensor_basic<T>();
631  for (size_type loc_jdof_bkp1 = 0, loc_ndof_bkp1 = grad_bkp1.size(); loc_jdof_bkp1 < loc_ndof_bkp1; ++loc_jdof_bkp1) {
632  for (size_type alpha = 0; alpha < d; ++alpha) {
633  size_type loc_jdof_bkp1d = d*loc_jdof_bkp1 + alpha;
634  for (size_type beta = 0; beta < d; ++beta) {
635  value[loc_idof](alpha,beta) += bar_a(loc_idof,loc_jdof_bkp1d)*grad_bkp1[loc_jdof_bkp1][beta];
636  }
637  }
638  }
639  }
640 }
641 // ----------------------------------------------------------------------------
642 // dofs for a scalar-valued function
643 // ----------------------------------------------------------------------------
644 // note: as virtual and template members are not available,
645 // the function "f" has been already evaluated on the hat_node[] set
646 // note2: moments on sides are scalar products with normals (not integrals on sides)
647 // for RT0 on the triangle we have the basis psi associated with momentum as:
648 // psi: matrix([x,y-1],[sqrt(2)*x,sqrt(2)*y],[x-1,y]);
649 // n : matrix([0,-1],[1/sqrt(2),1/sqrt(2)],[-1,0]);
650 // xm : matrix([1/2,0],[1/2,1/2],[0,1/2]);
651 // moment(i,f) := subst(xm[i][1],x,subst(xm[i][2],y,f.n[i]));
652 //
653 template<class T>
654 void
656  reference_element hat_K,
657  const Eigen::Matrix<value_type,Eigen::Dynamic,1>& f_xnod,
658  Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const
659 {
660 trace_macro ("_compute_dofs (hat_K="<<hat_K.name()<<")...");
661  base::_initialize_data_guard (hat_K);
662  size_type k = degree()-1;
663  size_type d = hat_K.dimension();
664  size_type loc_ndof = base::ndof (hat_K);
665 
666  dof.resize (loc_ndof);
667  if (d == 0) return;
668 
669  // side dofs are Lagrange ones, associated to nodes
670  size_type loc_inod = 0;
671  size_type loc_idof = 0;
672  for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(d-1); loc_isid < loc_nsid; ++loc_isid) {
673  reference_element hat_S = hat_K.subgeo (d-1, loc_isid);
674  point_basic<Float> hat_n;
675  Float hat_S_meas = hat_K.side_measure (loc_isid);
676  hat_K.side_normal (loc_isid, hat_n);
677  size_type loc_ndof_sid = _nnod_on_subgeo_remainder [d][hat_S.variant()];
678  for (size_type loc_idof_sid = 0; loc_idof_sid < loc_ndof_sid; ++loc_idof_sid) {
679  dof[loc_idof] = hat_S_meas*dot(f_xnod[loc_inod],hat_n);
680  loc_idof++;
681  loc_inod++;
682  }
683  }
684  // internal dofs
685  if (k == 0) return; // no internal dofs when k==0
686  const std::array<
687  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
688  ,3>&
689  bkm1_node_internal_d = _bkm1_node_internal_d [hat_K.variant()];
690 
691  size_type loc_ndof_boundary = loc_idof;
692  size_type loc_ndof_int_d = d*bkm1_node_internal_d[0].cols();
693  size_type first_loc_inod_int = loc_inod;
694  check_macro (loc_ndof == loc_ndof_boundary + loc_ndof_int_d,
695  "invalid internal dof count: loc_ndof="<<loc_ndof
696  << ", loc_ndof_boundary="<<loc_ndof_boundary
697  << ", loc_ndof_int_d="<<loc_ndof_int_d);
698 
699  size_type variant = hat_K.variant();
700  if (variant == reference_element::t) {
701  // interpolate all vector components on a lattice of internal nodes
702  for (size_type loc_idof_int = 0, loc_ndof_int = bkm1_node_internal_d[0].cols(); loc_idof_int < loc_ndof_int; ++loc_idof_int) {
703  size_type loc_inod = first_loc_inod_int + loc_idof_int;
704  for (size_type alpha = 0; alpha < d; ++alpha) {
705  dof [loc_idof] = f_xnod [loc_inod][alpha];
706  loc_idof++;
707  }
708  }
709  } else { // variant == qTPH : quadrature
710  for (size_type loc_idof_int = 0, loc_ndof_int = bkm1_node_internal_d[0].cols(); loc_idof_int < loc_ndof_int; ++loc_idof_int) {
711  for (size_type alpha = 0; alpha < d; ++alpha) {
712  loc_inod = first_loc_inod_int;
713  T sum = 0;
714  size_type inod_q = 0;
715  for (typename quadrature<T>::const_iterator iter_q = _quad.begin(hat_K),
716  last_q = _quad.end(hat_K); iter_q != last_q; iter_q++, inod_q++, ++loc_inod) {
717  sum += f_xnod [loc_inod][alpha]
718  *bkm1_node_internal_d[alpha] (inod_q, loc_idof_int)
719  *(*iter_q).w;
720  }
721  check_macro (loc_idof < loc_ndof, "invalid size");
722  dof [loc_idof] = sum;
723  loc_idof++;
724  }
725  }
726  }
727  check_macro (loc_idof == loc_ndof, "invalid dof count");
728 trace_macro ("_compute_dofs (hat_K="<<hat_K.name()<<") done");
729 }
730 // ----------------------------------------------------------------------------
731 // instanciation in library
732 // ----------------------------------------------------------------------------
733 #define _RHEOLEF_instanciation(T) \
734 template class basis_fem_RTk<T>;
735 
737 
738 }// namespace rheolef
rheolef::space_numbering::ndof
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:28
rheolef::reference_element::side_measure
Float side_measure(size_type loc_isid) const
Definition: reference_element.cc:213
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
warburton.icc
rheolef::reference_element::H
static const variant_type H
Definition: reference_element.h:81
rheolef::reference_element::n_node
static size_type n_node(variant_type variant, size_type order)
Definition: reference_element.cc:201
eigen_util.h
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
rheolef::point_basic
Definition: point.h:87
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::pointset_lagrange_equispaced
void pointset_lagrange_equispaced(reference_element hat_K, size_t order_in, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, size_t internal=0)
Definition: equispaced.icc:44
rheolef::basis_option::set_raw_polynomial
void set_raw_polynomial(raw_polynomial_type type)
Definition: basis_option.h:263
rheolef::dot
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
rheolef::reference_element::subgeo_local_vertex
size_type subgeo_local_vertex(size_type subgeo_dim, size_type loc_isid, size_type loc_jsidvert) const
Definition: reference_element.h:122
rheolef::details::basis_on_pointset_evaluate
void basis_on_pointset_evaluate(const Basis &b, const reference_element &hat_K, const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_x, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &vdm)
Definition: basis_on_pointset_evaluate.icc:31
rheolef::basis_fem_RTk::_initialize_cstor_sizes
void _initialize_cstor_sizes() const
Definition: basis_fem_RTk.cc:92
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::cond
T cond(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a)
Definition: eigen_util.h:62
rheolef::basis_raw_basic
Definition: basis_raw.h:107
rheolef::reference_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition: reference_element.h:110
rheolef::quadrature::const_iterator
rep::const_iterator const_iterator
Definition: quadrature.h:195
rheolef::value
rheolef::std value
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
rheolef::tensor_basic
Definition: tensor.h:90
rheolef::basis_option::warburton
Definition: basis_option.h:103
rheolef::basis_option::bernstein
Definition: basis_option.h:110
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::basis_option
see the basis_option page for the full documentation
Definition: basis_option.h:93
rheolef::basis_rep::size_type
reference_element::size_type size_type
Definition: basis.h:214
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::reference_element::subgeo
reference_element subgeo(size_type subgeo_dim, size_type loc_isid) const
Definition: reference_element.h:111
rheolef::basis_fem_RTk::grad_evaluate
void grad_evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &value) const
Definition: basis_fem_RTk.cc:603
rheolef::basis_fem_RTk::evaluate
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< value_type, Eigen::Dynamic, 1 > &value) const
Definition: basis_fem_RTk.cc:566
piola_fem_hdiv.h
mkgeo_ball.variant
variant
Definition: mkgeo_ball.sh:149
rheolef::basis_rep::_piola_fem
piola_fem< T > _piola_fem
Definition: basis.h:394
rheolef::basis_fem_RTk::hat_node
const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > & hat_node(reference_element hat_K) const
Definition: basis_fem_RTk.cc:71
a
Definition: diffusion_isotropic.h:25
rheolef::reference_element::name
char name() const
Definition: reference_element.h:100
rheolef::reference_element::variant
variant_type variant() const
Definition: reference_element.h:99
rheolef::basis_option::monomial
Definition: basis_option.h:109
rheolef::pointset_lagrange_warburton
void pointset_lagrange_warburton(reference_element hat_K, size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, bool map_on_reference_element=true)
Definition: warburton.icc:574
rheolef::basis_fem_RTk::_initialize_data
void _initialize_data(reference_element hat_K) const
Definition: basis_fem_RTk.cc:144
basis_fem_RTk.h
rheolef::basis_rep::_sopt
basis_option _sopt
Definition: basis.h:393
rheolef::basis_fem_RTk::inv_vdm
const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > & inv_vdm(reference_element hat_K) const
Definition: basis_fem_RTk.cc:85
rheolef::basis_rep::_name
std::string _name
Definition: basis.h:392
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
u
Definition: leveque.h:25
Float
see the Float page for the full documentation
basis_on_pointset_evaluate.icc
point
see the point page for the full documentation
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::basis_option::equispaced
Definition: basis_option.h:102
u
Float u(const point &x)
Definition: transmission_error.cc:26
rheolef::reference_element::dimension
size_type dimension() const
Definition: reference_element.h:101
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
rheolef::integrate_option::set_family
void set_family(family_type type)
Definition: integrate_option.h:260
mkgeo_ball.a
a
Definition: mkgeo_ball.sh:151
equispaced.icc
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::P
static const variant_type P
Definition: reference_element.h:80
rheolef::basis_fem_RTk::_compute_dofs
void _compute_dofs(reference_element hat_K, const Eigen::Matrix< value_type, Eigen::Dynamic, 1 > &f_xnod, Eigen::Matrix< T, Eigen::Dynamic, 1 > &dof) const
Definition: basis_fem_RTk.cc:655
rheolef::basis_fem_RTk::vdm
const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > & vdm(reference_element hat_K) const
Definition: basis_fem_RTk.cc:78
rheolef::basis_rep::standard_naming
static std::string standard_naming(std::string family_name, size_t degree, const basis_option &sopt)
Definition: basis_rep.cc:44
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::basis_fem_RTk::basis_fem_RTk
basis_fem_RTk(size_type k, const basis_option &sopt)
Definition: basis_fem_RTk.cc:43
rheolef::invert
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition: tiny_lu.h:127
rheolef::piola_fem_hdiv
Definition: piola_fem_hdiv.h:74
rheolef::reference_element::side_normal
void side_normal(size_type loc_isid, point_basic< Float > &hat_n) const
Definition: reference_element.cc:232
rheolef::basis_option::dubiner
Definition: basis_option.h:111
epsilon
Float epsilon
Definition: transmission_error.cc:25
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
rheolef::reference_element::vertex
const point_basic< Float > & vertex(size_type iloc) const
Definition: reference_element.cc:251
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
rheolef::quadrature
Definition: quadrature.h:185
rheolef::reference_element::n_vertex
size_type n_vertex() const
Definition: reference_element.h:103
rheolef::basis_fem_RTk::~basis_fem_RTk
~basis_fem_RTk()
Definition: basis_fem_RTk.cc:39
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
basis_ordering.icc
rk::beta
Float beta[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:60
rheolef::basis_fem_RTk::_b_pre_kp1
basis_raw_basic< T > _b_pre_kp1
Definition: basis_fem_RTk.h:120
rheolef::basis_fem_RTk::family_name
std::string family_name() const
Definition: basis_fem_RTk.h:71
rheolef::basis_rep
Definition: basis.h:209
T
Expr1::float_type T
Definition: field_expr.h:218
trace_macro
#define trace_macro(message)
Definition: dis_macros.h:111