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