GetFEM  5.4.2
bgeot_geometric_trans.cc
1 /*===========================================================================
2 
3  Copyright (C) 2000-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 
23 #include "getfem/dal_singleton.h"
24 #include "getfem/dal_tree_sorted.h"
27 #include "getfem/bgeot_torus.h"
28 
29 namespace bgeot {
30 
31  std::vector<scalar_type>& __aux1(){
32  THREAD_SAFE_STATIC std::vector<scalar_type> v;
33  return v;
34  }
35 
36  std::vector<scalar_type>& __aux2(){
37  THREAD_SAFE_STATIC std::vector<scalar_type> v;
38  return v;
39  }
40 
41  std::vector<scalar_type>& __aux3(){
42  THREAD_SAFE_STATIC std::vector<scalar_type> v;
43  return v;
44  }
45 
46  std::vector<long>& __ipvt_aux(){
47  THREAD_SAFE_STATIC std::vector<long> vi;
48  return vi;
49  }
50 
51  // Optimized matrix mult for small matrices. To be verified.
52  // Multiply the matrix A of size MxN by B of size NxP in C of size MxP
53  void mat_mult(const scalar_type *A, const scalar_type *B, scalar_type *C,
54  size_type M, size_type N, size_type P) {
55  if (N != 0) {
56  auto itC = C; auto itB = B;
57  for (size_type j = 0; j < P; ++j, itB += N)
58  for (size_type i = 0; i < M; ++i, ++itC) {
59  auto itA = A+i, itB1 = itB;
60  *itC = (*itA) * (*itB1);
61  for (size_type k = 1; k < N; ++k)
62  { itA += M; ++itB1; *itC += (*itA) * (*itB1); }
63  }
64  } else std::fill(C, C+M*P, scalar_type(0));
65  }
66 
67  // Optimized matrix mult for small matrices.
68  // Multiply the matrix A of size MxN by the transpose of B of size PxN
69  // in C of size MxP
70  void mat_tmult(const scalar_type *A, const scalar_type *B, scalar_type *C,
71  size_type M, size_type N, size_type P) {
72  auto itC = C;
73  switch (N) {
74  case 0 : std::fill(C, C+M*P, scalar_type(0)); break;
75  case 1 :
76  for (size_type j = 0; j < P; ++j)
77  for (size_type i = 0; i < M; ++i, ++itC) {
78  auto itA = A+i, itB = B+j;
79  *itC = (*itA) * (*itB);
80  }
81  break;
82  case 2 :
83  for (size_type j = 0; j < P; ++j)
84  for (size_type i = 0; i < M; ++i, ++itC) {
85  auto itA = A+i, itB = B+j;
86  *itC = (*itA) * (*itB);
87  itA += M; itB += P; *itC += (*itA) * (*itB);
88  }
89  break;
90  case 3 :
91  for (size_type j = 0; j < P; ++j)
92  for (size_type i = 0; i < M; ++i, ++itC) {
93  auto itA = A+i, itB = B+j;
94  *itC = (*itA) * (*itB);
95  itA += M; itB += P; *itC += (*itA) * (*itB);
96  itA += M; itB += P; *itC += (*itA) * (*itB);
97  }
98  break;
99  default :
100  for (size_type j = 0; j < P; ++j)
101  for (size_type i = 0; i < M; ++i, ++itC) {
102  auto itA = A+i, itB = B+j;
103  *itC = (*itA) * (*itB);
104  for (size_type k = 1; k < N; ++k)
105  { itA += M; itB += P; *itC += (*itA) * (*itB); }
106  }
107  }
108  }
109 
110 
111 
112 
113  // Optimized lu_factor for small square matrices
114  size_type lu_factor(scalar_type *A, std::vector<long> &ipvt,
115  size_type N) {
116  size_type info(0), i, j, jp, N_1 = N-1;
117 
118  if (N) {
119  for (j = 0; j < N_1; ++j) {
120  auto it = A + (j*(N+1));
121  scalar_type max = gmm::abs(*it); jp = j;
122  for (i = j+1; i < N; ++i) {
123  scalar_type ap = gmm::abs(*(++it));
124  if (ap > max) { jp = i; max = ap; }
125  }
126  ipvt[j] = long(jp + 1);
127 
128  if (max == scalar_type(0)) { info = j + 1; break; }
129  if (jp != j) {
130  auto it1 = A+jp, it2 = A+j;
131  for (i = 0; i < N; ++i, it1+=N, it2+=N) std::swap(*it1, *it2);
132  }
133  it = A + (j*(N+1)); max = *it++;
134  for (i = j+1; i < N; ++i) *it++ /= max;
135  auto it22 = A + (j*N + j+1), it11 = it22;
136  auto it3 = A + ((j+1)*N+j);
137  for (size_type l = j+1; l < N; ++l) {
138  it11 += N;
139  auto it1 = it11, it2 = it22;
140  scalar_type a = *it3; it3 += N;
141  for (size_type k = j+1; k < N; ++k) *it1++ -= *it2++ * a;
142  }
143 
144  }
145  ipvt[N-1] = long(N);
146  }
147  return info;
148  }
149 
150  static void lower_tri_solve(const scalar_type *T, scalar_type *x, int N,
151  bool is_unit) {
152  scalar_type x_j;
153  for (int j = 0; j < N; ++j) {
154  auto itc = T + j*N, it = itc+(j+1), ite = itc+N;
155  auto itx = x + j;
156  if (!is_unit) *itx /= itc[j];
157  x_j = *itx++;
158  for (; it != ite ; ++it, ++itx) *itx -= x_j * (*it);
159  }
160  }
161 
162  static void upper_tri_solve(const scalar_type *T, scalar_type *x, int N,
163  bool is_unit) {
164  scalar_type x_j;
165  for (int j = N - 1; j >= 0; --j) {
166  auto itc = T + j*N, it = itc, ite = itc+j;
167  auto itx = x;
168  if (!is_unit) x[j] /= itc[j];
169  for (x_j = x[j]; it != ite ; ++it, ++itx) *itx -= x_j * (*it);
170  }
171  }
172 
173  static void lu_solve(const scalar_type *LU, const std::vector<long> &ipvt,
174  scalar_type *x, scalar_type *b, int N) {
175  std::copy(b, b+N, x);
176  for(int i = 0; i < N; ++i)
177  { long perm = ipvt[i]-1; if(i != perm) std::swap(x[i], x[perm]); }
178  bgeot::lower_tri_solve(LU, x, N, true);
179  bgeot::upper_tri_solve(LU, x, N, false);
180  }
181 
182  scalar_type lu_det(const scalar_type *LU, const std::vector<long> &ipvt,
183  size_type N) {
184  scalar_type det(1);
185  for (size_type j = 0; j < N; ++j) det *= *(LU+j*(N+1));
186  for(long i = 0; i < long(N); ++i) if (i != ipvt[i]-1) { det = -det; }
187  return det;
188  }
189 
190  scalar_type lu_det(const scalar_type *A, size_type N) {
191  switch (N) {
192  case 1: return *A;
193  case 2: return (*A) * (A[3]) - (A[1]) * (A[2]);
194  case 3:
195  {
196  scalar_type a0 = A[4]*A[8] - A[5]*A[7], a3 = A[5]*A[6] - A[3]*A[8];
197  scalar_type a6 = A[3]*A[7] - A[4]*A[6];
198  return A[0] * a0 + A[1] * a3 + A[2] * a6;
199  }
200  default:
201  {
202  size_type NN = N*N;
203  if (__aux1().size() < NN) __aux1().resize(N*N);
204  std::copy(A, A+NN, __aux1().begin());
205  __ipvt_aux().resize(N);
206  lu_factor(&(*(__aux1().begin())), __ipvt_aux(), N);
207  return lu_det(&(*(__aux1().begin())), __ipvt_aux(), N);
208  }
209  }
210  }
211 
212  void lu_inverse(const scalar_type *LU, const std::vector<long> &ipvt,
213  scalar_type *A, size_type N) {
214  __aux2().resize(N); gmm::clear(__aux2());
215  __aux3().resize(N);
216  for(size_type i = 0; i < N; ++i) {
217  __aux2()[i] = scalar_type(1);
218  bgeot::lu_solve(LU, ipvt, A+i*N, &(*(__aux2().begin())), int(N));
219  __aux2()[i] = scalar_type(0);
220  }
221  }
222 
223  scalar_type lu_inverse(scalar_type *A, size_type N, bool doassert) {
224  switch (N) {
225  case 1:
226  {
227  scalar_type det = *A;
228  GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix");
229  *A = scalar_type(1)/det;
230  return det;
231  }
232  case 2:
233  {
234  scalar_type a = *A, b = A[2], c = A[1], d = A[3];
235  scalar_type det = a * d - b * c;
236  GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix");
237  *A++ = d/det; *A++ /= -det; *A++ /= -det; *A = a/det;
238  return det;
239  }
240  case 3:
241  {
242  scalar_type a0 = A[4]*A[8] - A[5]*A[7], a1 = A[5]*A[6] - A[3]*A[8];
243  scalar_type a2 = A[3]*A[7] - A[4]*A[6];
244  scalar_type det = A[0] * a0 + A[1] * a1 + A[2] * a2;
245  GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix");
246  scalar_type a3 = (A[2]*A[7] - A[1]*A[8]), a6 = (A[1]*A[5] - A[2]*A[4]);
247  scalar_type a4 = (A[0]*A[8] - A[2]*A[6]), a7 = (A[2]*A[3] - A[0]*A[5]);
248  scalar_type a5 = (A[1]*A[6] - A[0]*A[7]), a8 = (A[0]*A[4] - A[1]*A[3]);
249  *A++ = a0 / det; *A++ = a3 / det; *A++ = a6 / det;
250  *A++ = a1 / det; *A++ = a4 / det; *A++ = a7 / det;
251  *A++ = a2 / det; *A++ = a5 / det; *A++ = a8 / det;
252  return det;
253  }
254  default:
255  {
256  size_type NN = N*N;
257  if (__aux1().size() < NN) __aux1().resize(NN);
258  std::copy(A, A+NN, __aux1().begin());
259  __ipvt_aux().resize(N);
260  size_type info = lu_factor(&(*(__aux1().begin())), __ipvt_aux(), N);
261  if (doassert) GMM_ASSERT1(!info, "Non invertible matrix, pivot = "
262  << info);
263  if (!info) lu_inverse(&(*(__aux1().begin())), __ipvt_aux(), A, N);
264  return lu_det(&(*(__aux1().begin())), __ipvt_aux(), N);
265  }
266  }
267  }
268 
269 
270 
272  (const base_matrix &G, const base_matrix &pc, base_matrix &K) const {
273  // gmm::mult(G, pc, K);
274  // Faster than the lapack call on my config ...
275  size_type N=gmm::mat_nrows(G), P=gmm::mat_nrows(pc), Q=gmm::mat_ncols(pc);
276  if (N && P && Q) {
277  auto itK = K.begin();
278  for (size_type j = 0; j < Q; ++j) {
279  auto itpc_j = pc.begin() + j*P, itG_b = G.begin();
280  for (size_type i = 0; i < N; ++i, ++itG_b) {
281  auto itG = itG_b, itpc = itpc_j;
282  scalar_type a = *(itG) * (*itpc);
283  for (size_type k = 1; k < P; ++k)
284  { itG += N; a += *(itG) * (*++itpc); }
285  *itK++ = a;
286  }
287  }
288  } else gmm::clear(K);
289  }
290 
292  if (!have_xref()) {
293  if (pspt_) xref_ = (*pspt_)[ii_];
294  else GMM_ASSERT1(false, "Missing xref");
295  }
296  return xref_;
297  }
298 
300  if (!have_xreal()) {
301  if (have_pgp()) {
302  xreal_ = pgp_->transform(ii_, G());
303  } else xreal_ = pgt()->transform(xref(),G());
304  }
305  return xreal_;
306  }
307 
309  GMM_ASSERT1(have_G(),
310  "Convex center can be provided only if matrix G is available");
311  if (!have_cv_center_) {
312  cv_center_.resize(G().nrows());
313  size_type nb_pts = G().ncols();
314  for (size_type i=0; i < nb_pts; i++)
315  gmm::add(gmm::mat_col(G(),i), cv_center_);
316  gmm::scale(cv_center_, scalar_type(1)/scalar_type(nb_pts));
317  have_cv_center_ = true;
318  }
319  return cv_center_;
320  }
321 
322  void geotrans_interpolation_context::compute_J() const {
323  GMM_ASSERT1(have_G() && have_pgt(), "Unable to compute J\n");
324  size_type P = pgt_->structure()->dim();
325  const base_matrix &KK = K();
326  if (P != N()) {
327  B_factors.base_resize(P, P);
328  gmm::mult(gmm::transposed(KK), KK, B_factors);
329  // gmm::abs below because on flat convexes determinant could be -1e-27.
330  J__ = J_ =::sqrt(gmm::abs(bgeot::lu_inverse(&(*(B_factors.begin())), P)));
331  } else {
332  auto it = &(*(KK.begin()));
333  switch (P) {
334  case 1: J__ = *it; break;
335  case 2: J__ = (*it) * (it[3]) - (it[1]) * (it[2]); break;
336  case 3:
337  {
338  B_.base_resize(P, P); // co-factors
339  auto itB = B_.begin();
340  scalar_type a0 = itB[0] = it[4]*it[8] - it[5]*it[7];
341  scalar_type a1 = itB[1] = it[5]*it[6] - it[3]*it[8];
342  scalar_type a2 = itB[2] = it[3]*it[7] - it[4]*it[6];
343  J__ = it[0] * a0 + it[1] * a1 + it[2] * a2;
344  } break;
345  default:
346  B_factors.base_resize(P, P); // store factorization for B computation
347  gmm::copy(gmm::transposed(KK), B_factors);
348  ipvt.resize(P);
349  bgeot::lu_factor(&(*(B_factors.begin())), ipvt, P);
350  J__ = bgeot::lu_det(&(*(B_factors.begin())), ipvt, P);
351  break;
352  }
353  J_ = gmm::abs(J__);
354  }
355  have_J_ = true;
356  }
357 
358  const base_matrix& geotrans_interpolation_context::K() const {
359  if (!have_K()) {
360  GMM_ASSERT1(have_G() && have_pgt(), "Unable to compute K\n");
361  size_type P = pgt_->structure()->dim();
362  K_.base_resize(N(), P);
363  if (have_pgp()) {
364  pgt_->compute_K_matrix(*G_, pgp_->grad(ii_), K_);
365  } else {
366  PC.base_resize(pgt_->nb_points(), P);
367  pgt_->poly_vector_grad(xref(), PC);
368  pgt_->compute_K_matrix(*G_, PC, K_);
369  }
370  have_K_ = true;
371  }
372  return K_;
373  }
374 
375  const base_matrix& geotrans_interpolation_context::B() const {
376  if (!have_B()) {
377  const base_matrix &KK = K();
378  size_type P = pgt_->structure()->dim(), N_ = gmm::mat_nrows(KK);
379  B_.base_resize(N_, P);
380  if (!have_J_) compute_J();
381  GMM_ASSERT1(J__ != scalar_type(0), "Non invertible matrix");
382  if (P != N_) {
383  gmm::mult(KK, B_factors, B_);
384  } else {
385  switch(P) {
386  case 1: B_(0, 0) = scalar_type(1) / J__; break;
387  case 2:
388  {
389  auto it = &(*(KK.begin())); auto itB = &(*(B_.begin()));
390  *itB++ = it[3] / J__;
391  *itB++ = -it[2] / J__;
392  *itB++ = -it[1] / J__;
393  *itB = (*it) / J__;
394  } break;
395  case 3:
396  {
397  auto it = &(*(KK.begin())); auto itB = &(*(B_.begin()));
398  *itB++ /= J__; *itB++ /= J__; *itB++ /= J__;
399  *itB++ = (it[2]*it[7] - it[1]*it[8]) / J__;
400  *itB++ = (it[0]*it[8] - it[2]*it[6]) / J__;
401  *itB++ = (it[1]*it[6] - it[0]*it[7]) / J__;
402  *itB++ = (it[1]*it[5] - it[2]*it[4]) / J__;
403  *itB++ = (it[2]*it[3] - it[0]*it[5]) / J__;
404  *itB = (it[0]*it[4] - it[1]*it[3]) / J__;
405  } break;
406  default:
407  bgeot::lu_inverse(&(*(B_factors.begin())), ipvt, &(*(B_.begin())), P);
408  break;
409  }
410  }
411  have_B_ = true;
412  }
413  return B_;
414  }
415 
416  const base_matrix& geotrans_interpolation_context::B3() const {
417  if (!have_B3()) {
418  const base_matrix &BB = B();
419  size_type P=gmm::mat_ncols(BB), N_=gmm::mat_nrows(BB);
420  B3_.base_resize(N_*N_, P*P);
421  for (short_type i = 0; i < P; ++i)
422  for (short_type j = 0; j < P; ++j)
423  for (short_type k = 0; k < N_; ++k)
424  for (short_type l = 0; l < N_; ++l)
425  B3_(k + N_*l, i + P*j) = BB(k, i) * BB(l, j);
426  have_B3_ = true;
427  }
428  return B3_;
429  }
430 
431  const base_matrix& geotrans_interpolation_context::B32() const {
432  if (!have_B32()) {
433  const base_matrix &BB = B();
434  size_type P=gmm::mat_ncols(BB), N_=gmm::mat_nrows(BB);
435  B32_.base_resize(N_*N_, P);
436  if (!pgt()->is_linear()) {
437  base_matrix B2(P*P, P), Htau(N_, P*P);
438  if (have_pgp()) {
439  gmm::mult(G(), pgp_->hessian(ii_), Htau);
440  } else {
441  /* very inefficient of course... */
442  PC.base_resize(pgt()->nb_points(), P*P);
443  pgt()->poly_vector_hess(xref(), PC);
444  gmm::mult(G(), PC, Htau);
445  }
446  for (short_type i = 0; i < P; ++i)
447  for (short_type j = 0; j < P; ++j)
448  for (short_type k = 0; k < P; ++k)
449  for (short_type l = 0; l < N_; ++l)
450  B2(i + P*j, k) += Htau(l, i + P*j) * BB(l,k);
451  gmm::mult(B3(), B2, B32_);
452  } else gmm::clear(B32_);
453  have_B32_ = true;
454  }
455  return B32_;
456  }
457 
459  xref_ = P;
460  if (pgt_ && !pgt()->is_linear())
461  { have_K_ = have_B_ = have_B3_ = have_B32_ = have_J_ = false; }
462  xreal_.resize(0); ii_ = size_type(-1); pspt_ = 0;
463  }
464 
465 
467  const base_matrix &G) const {
468  size_type N = G.nrows(), k = nb_points();
469  base_node P(N); base_vector val(k);
470  poly_vector_val(pt, val);
471  base_matrix::const_iterator git = G.begin();
472  for (size_type l = 0; l < k; ++l) {
473  scalar_type a = val[l];
474  base_node::iterator pit = P.begin(), pite = P.end();
475  for (; pit != pite; ++git, ++pit) *pit += a * (*git);
476  }
477  return P;
478  }
479 
480  void geometric_trans::fill_standard_vertices() {
481  vertices_.resize(0);
482  for (size_type ip = 0; ip < nb_points(); ++ip) {
483  bool vertex = true;
484  for (size_type i = 0; i < cvr->points()[ip].size(); ++i)
485  if (gmm::abs(cvr->points()[ip][i]) > 1e-10
486  && gmm::abs(cvr->points()[ip][i]-1.0) > 1e-10
487  && gmm::abs(cvr->points()[ip][i]+1.0) > 1e-10)
488  { vertex = false; break; }
489  if (vertex) vertices_.push_back(ip);
490  }
491  dim_type dimension = dim();
492  if (dynamic_cast<const torus_geom_trans *>(this)) --dimension;
493  GMM_ASSERT1(vertices_.size() > dimension, "Internal error");
494  }
495 
496  /* ******************************************************************** */
497  /* Instantied geometric transformations. */
498  /* ******************************************************************** */
499 
500  template <class FUNC>
501  struct igeometric_trans : public geometric_trans {
502 
503  std::vector<FUNC> trans;
504  mutable std::vector<std::vector<FUNC>> grad_, hess_;
505  mutable bool grad_computed_ = false;
506  mutable bool hess_computed_ = false;
507 
508  void compute_grad_() const {
509  if (grad_computed_) return;
510  GLOBAL_OMP_GUARD
511  if (grad_computed_) return;
512  size_type R = trans.size();
513  dim_type n = dim();
514  grad_.resize(R);
515  for (size_type i = 0; i < R; ++i) {
516  grad_[i].resize(n);
517  for (dim_type j = 0; j < n; ++j) {
518  grad_[i][j] = trans[i]; grad_[i][j].derivative(j);
519  }
520  }
521  grad_computed_ = true;
522  }
523 
524  void compute_hess_() const {
525  if (hess_computed_) return;
526  GLOBAL_OMP_GUARD
527  if (hess_computed_) return;
528  size_type R = trans.size();
529  dim_type n = dim();
530  hess_.resize(R);
531  for (size_type i = 0; i < R; ++i) {
532  hess_[i].resize(n*n);
533  for (dim_type j = 0; j < n; ++j) {
534  for (dim_type k = 0; k < n; ++k) {
535  hess_[i][j+k*n] = trans[i];
536  hess_[i][j+k*n].derivative(j); hess_[i][j+k*n].derivative(k);
537  }
538  }
539  }
540  hess_computed_ = true;
541  }
542 
543  virtual void poly_vector_val(const base_node &pt, base_vector &val) const {
544  val.resize(nb_points());
545  for (size_type k = 0; k < nb_points(); ++k)
546  val[k] = to_scalar(trans[k].eval(pt.begin()));
547  }
548 
549  virtual void poly_vector_val(const base_node &pt,
550  const convex_ind_ct &ind_ct,
551  base_vector &val) const {
552  size_type nb_funcs=ind_ct.size();
553  val.resize(nb_funcs);
554  for (size_type k = 0; k < nb_funcs; ++k)
555  val[k] = to_scalar(trans[ind_ct[k]].eval(pt.begin()));
556  }
557 
558  virtual void poly_vector_grad(const base_node &pt, base_matrix &pc) const {
559  if (!grad_computed_) compute_grad_();
560  FUNC PP;
561  pc.base_resize(nb_points(),dim());
562  for (size_type i = 0; i < nb_points(); ++i)
563  for (dim_type n = 0; n < dim(); ++n)
564  pc(i, n) = to_scalar(grad_[i][n].eval(pt.begin()));
565  }
566 
567  virtual void poly_vector_grad(const base_node &pt,
568  const convex_ind_ct &ind_ct,
569  base_matrix &pc) const {
570  if (!grad_computed_) compute_grad_();
571  FUNC PP;
572  size_type nb_funcs=ind_ct.size();
573  pc.base_resize(nb_funcs,dim());
574  for (size_type i = 0; i < nb_funcs; ++i)
575  for (dim_type n = 0; n < dim(); ++n)
576  pc(i, n) = to_scalar(grad_[ind_ct[i]][n].eval(pt.begin()));
577  }
578 
579  virtual void poly_vector_hess(const base_node &pt, base_matrix &pc) const {
580  if (!hess_computed_) compute_hess_();
581  FUNC PP, QP;
582  pc.base_resize(nb_points(),dim()*dim());
583  for (size_type i = 0; i < nb_points(); ++i)
584  for (dim_type n = 0; n < dim(); ++n) {
585  for (dim_type m = 0; m <= n; ++m)
586  pc(i, n*dim()+m) = pc(i, m*dim()+n) =
587  to_scalar(hess_[i][m*dim()+n].eval(pt.begin()));
588  }
589  }
590 
591  };
592 
593  typedef igeometric_trans<base_poly> poly_geometric_trans;
594  typedef igeometric_trans<polynomial_composite> comppoly_geometric_trans;
595  typedef igeometric_trans<base_rational_fraction> fraction_geometric_trans;
596 
597  /* ******************************************************************** */
598  /* transformation on simplex. */
599  /* ******************************************************************** */
600 
601  struct simplex_trans_ : public poly_geometric_trans {
602  void calc_base_func(base_poly &p, size_type i, short_type K) const {
603  dim_type N = dim();
604  base_poly l0(N, 0), l1(N, 0);
605  power_index w(short_type(N+1));
606  l0.one(); l1.one(); p = l0;
607  for (short_type nn = 0; nn < N; ++nn) l0 -= base_poly(N, 1, nn);
608 
609  w[0] = K;
610  for (int nn = 1; nn <= N; ++nn) {
611  w[nn]=short_type(floor(0.5+(((cvr->points())[i])[nn-1]*double(K))));
612  w[0]=short_type(w[0]-w[nn]);
613  }
614 
615  for (short_type nn = 0; nn <= N; ++nn)
616  for (short_type j = 0; j < w[nn]; ++j)
617  if (nn == 0)
618  p *= (l0 * (scalar_type(K) / scalar_type(j+1)))
619  - (l1 * (scalar_type(j) / scalar_type(j+1)));
620  else
621  p *= (base_poly(N, 1, short_type(nn-1)) * (scalar_type(K) / scalar_type(j+1)))
622  - (l1 * (scalar_type(j) / scalar_type(j+1)));
623  }
624 
625  simplex_trans_(dim_type nc, short_type k) {
626  cvr = simplex_of_reference(nc, k);
627  size_type R = cvr->structure()->nb_points();
628  is_lin = (k == 1);
629  complexity_ = k;
630  trans.resize(R);
631  for (size_type r = 0; r < R; ++r) calc_base_func(trans[r], r, k);
632  fill_standard_vertices();
633  }
634  };
635 
636  static pgeometric_trans
637  PK_gt(gt_param_list &params,
638  std::vector<dal::pstatic_stored_object> &dependencies) {
639  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
640  << params.size() << " should be 2.");
641  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
642  "Bad type of parameters");
643  int n = int(::floor(params[0].num() + 0.01));
644  int k = int(::floor(params[1].num() + 0.01));
645  GMM_ASSERT1(n >= 0 && n < 100 && k >= 0 && k <= 150 &&
646  double(n) == params[0].num() && double(k) == params[1].num(),
647  "Bad parameters");
648  dependencies.push_back(simplex_of_reference(dim_type(n), dim_type(k)));
649  return std::make_shared<simplex_trans_>(dim_type(n), dim_type(k));
650  }
651 
652  /* ******************************************************************** */
653  /* direct product transformation */
654  /* ******************************************************************** */
655 
656  struct cv_pr_t_ : public poly_geometric_trans {
657  cv_pr_t_(const poly_geometric_trans *a, const poly_geometric_trans *b) {
658  cvr = convex_ref_product(a->convex_ref(), b->convex_ref());
659  is_lin = false;
660  complexity_ = a->complexity() * b->complexity();
661 
662  size_type n1 = a->nb_points(), n2 = b->nb_points();
663  trans.resize(n1 * n2);
664  for (size_type i1 = 0; i1 < n1; ++i1)
665  for (size_type i2 = 0; i2 < n2; ++i2) {
666  trans[i1 + i2 * n1] = a->trans[i1];
667  trans[i1 + i2 * n1].direct_product(b->trans[i2]);
668  }
669  for (size_type i2 = 0; i2 < b->nb_vertices(); ++i2)
670  for (size_type i1 = 0; i1 < a->nb_vertices(); ++i1)
671  vertices_.push_back(a->vertices()[i1] + b->vertices()[i2] * n1);
672  }
673  };
674 
675  static pgeometric_trans product_gt(gt_param_list &params,
676  std::vector<dal::pstatic_stored_object> &dependencies) {
677  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
678  << params.size() << " should be 2.");
679  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
680  "Bad type of parameters");
681  pgeometric_trans a = params[0].method();
682  pgeometric_trans b = params[1].method();
683  dependencies.push_back(a); dependencies.push_back(b);
684  dependencies.push_back(convex_ref_product(a->convex_ref(),
685  b->convex_ref()));
686  const poly_geometric_trans *aa
687  = dynamic_cast<const poly_geometric_trans *>(a.get());
688  const poly_geometric_trans *bb
689  = dynamic_cast<const poly_geometric_trans *>(b.get());
690  GMM_ASSERT1(aa && bb, "The product of geometric transformations "
691  "is only defined for polynomial ones");
692  return std::make_shared<cv_pr_t_>(aa, bb);
693  }
694 
695  /* ******************************************************************** */
696  /* linear direct product transformation. */
697  /* ******************************************************************** */
698 
699  struct cv_pr_tl_ : public poly_geometric_trans {
700  cv_pr_tl_(const poly_geometric_trans *a, const poly_geometric_trans *b) {
701  GMM_ASSERT1(a->is_linear() && b->is_linear(),
702  "linear product of non-linear transformations");
703  cvr = convex_ref_product(a->convex_ref(), b->convex_ref());
704  is_lin = true;
705  complexity_ = std::max(a->complexity(), b->complexity());
706 
707  trans.resize(a->nb_points() * b->nb_points());
708  std::fill(trans.begin(), trans.end(), null_poly(dim()));
709 
710  std::stringstream name;
711  name << "GT_PK(" << int(dim()) << ",1)";
712  pgeometric_trans pgt_ = geometric_trans_descriptor(name.str());
713  const poly_geometric_trans *pgt
714  = dynamic_cast<const poly_geometric_trans *>(pgt_.get());
715 
716  for (size_type i = 0; i <= dim(); ++i)
717  trans[cvr->structure()->ind_dir_points()[i]]
718  = pgt->trans[i];
719  for (size_type i2 = 0; i2 < b->nb_vertices(); ++i2)
720  for (size_type i1 = 0; i1 < a->nb_vertices(); ++i1)
721  vertices_.push_back(a->vertices()[i1]
722  + b->vertices()[i2] * a->nb_points());
723  }
724  };
725 
726  static pgeometric_trans linear_product_gt(gt_param_list &params,
727  std::vector<dal::pstatic_stored_object> &dependencies) {
728  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
729  << params.size() << " should be 2.");
730  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
731  "Bad type of parameters");
732  pgeometric_trans a = params[0].method();
733  pgeometric_trans b = params[1].method();
734  dependencies.push_back(a); dependencies.push_back(b);
735  dependencies.push_back(convex_ref_product(a->convex_ref(),
736  b->convex_ref()));
737  const poly_geometric_trans *aa
738  = dynamic_cast<const poly_geometric_trans *>(a.get());
739  const poly_geometric_trans *bb
740  = dynamic_cast<const poly_geometric_trans *>(b.get());
741  GMM_ASSERT1(aa && bb, "The product of geometric transformations "
742  "is only defined for polynomial ones");
743  return std::make_shared<cv_pr_tl_>(aa, bb);
744  }
745 
746  /* ******************************************************************** */
747  /* parallelepiped transformation. */
748  /* ******************************************************************** */
749 
750  static pgeometric_trans QK_gt(gt_param_list &params,
751  std::vector<dal::pstatic_stored_object> &) {
752  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
753  << params.size() << " should be 2.");
754  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
755  "Bad type of parameters");
756  int n = int(::floor(params[0].num() + 0.01));
757  int k = int(::floor(params[1].num() + 0.01));
758  GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
759  double(n) == params[0].num() && double(k) == params[1].num(),
760  "Bad parameters");
761  std::stringstream name;
762  if (n == 1)
763  name << "GT_PK(1," << k << ")";
764  else
765  name << "GT_PRODUCT(GT_QK(" << n-1 << "," << k << "),GT_PK(1,"
766  << k << "))";
767  return geometric_trans_descriptor(name.str());
768  }
769 
770  static pgeometric_trans
771  prism_pk_gt(gt_param_list &params,
772  std::vector<dal::pstatic_stored_object> &) {
773  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
774  << params.size() << " should be 2.");
775  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
776  "Bad type of parameters");
777  int n = int(::floor(params[0].num() + 0.01));
778  int k = int(::floor(params[1].num() + 0.01));
779  GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
780  double(n) == params[0].num() && double(k) == params[1].num(),
781  "Bad parameters");
782  std::stringstream name;
783  name << "GT_PRODUCT(GT_PK(" << n-1 << "," << k << "),GT_PK(1,"
784  << k << "))";
785  return geometric_trans_descriptor(name.str());
786  }
787 
788  static pgeometric_trans
789  linear_qk(gt_param_list &params,
790  std::vector<dal::pstatic_stored_object> &) {
791  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
792  << params.size() << " should be 1.");
793  GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
794  int n = int(::floor(params[0].num() + 0.01));
795  return parallelepiped_linear_geotrans(n);
796  }
797 
798 
799  /* ******************************************************************** */
800  /* Incomplete Q2 geometric transformation for n=2 or 3. */
801  /* ******************************************************************** */
802  /* By Yao Koutsawa <[email protected]> 2012-12-10 */
803 
804  struct Q2_incomplete_trans_: public poly_geometric_trans {
805  Q2_incomplete_trans_(dim_type nc) {
806  cvr = Q2_incomplete_of_reference(nc);
807  size_type R = cvr->structure()->nb_points();
808  is_lin = false;
809  complexity_ = 2;
810  trans.resize(R);
811 
812  if (nc == 2) {
813  std::stringstream s
814  ( "1 - 2*x^2*y - 2*x*y^2 + 2*x^2 + 5*x*y + 2*y^2 - 3*x - 3*y;"
815  "4*(x^2*y - x^2 - x*y + x);"
816  "2*x*y*y - 2*x*x*y + 2*x*x - x*y - x;"
817  "4*(x*y*y - x*y - y*y + y);"
818  "4*(x*y - x*y*y);"
819  "2*x*x*y - 2*x*y*y - x*y + 2*y*y - y;"
820  "4*(x*y - x*x*y);"
821  "2*x*x*y + 2*x*y*y - 3*x*y;");
822 
823  for (int i = 0; i < 8; ++i)
824  trans[i] = read_base_poly(2, s);
825  } else {
826  std::stringstream s
827  ("1 + 2*x^2*y*z + 2*x*y^2*z + 2*x*y*z^2"
828  " - 2*x^2*y - 2*x^2*z - 2*x*y^2 - 2*y^2*z - 2*y*z^2 - 2*x*z^2 - 7*x*y*z"
829  " + 2*x^2 + 2*y^2 + 2*z^2 + 5*y*z + 5*x*z + 5*x*y - 3*x - 3*y - 3*z;"
830  "4*( - x^2*y*z + x*y*z + x^2*z - x*z + x^2*y - x*y - x^2 + x);"
831  "2*x^2*y*z - 2*x*y^2*z - 2*x*y*z^2"
832  " - 2*x^2*y - 2*x^2*z + 2*x*y^2 + 2*x*z^2 + 3*x*y*z + 2*x^2 - x*y - x*z - x;"
833  "4*( - x*y^2*z + x*y^2 + y^2*z + x*y*z - x*y - y^2 - y*z + y);"
834  "4*(x*y^2*z - x*y^2 - x*y*z + x*y);"
835  " - 2*x^2*y*z + 2*x*y^2*z - 2*x*y*z^2"
836  " + 2*x^2*y - 2*x*y^2 - 2*y^2*z + 2*y*z^2 + 3*x*y*z - x*y + 2*y^2 - y*z - y;"
837  "4*(x^2*y*z - x^2*y - x*y*z + x*y);"
838  " - 2*x^2*y*z - 2*x*y^2*z + 2*x*y*z^2 + 2*x^2*y + 2*x*y^2 + x*y*z - 3*x*y;"
839  "4*( - x*y*z^2 + x*z^2 + y*z^2 + x*y*z - x*z - y*z - z^2 + z);"
840  "4*(x*y*z^2 - x*y*z - x*z^2 + x*z);"
841  "4*(x*y*z^2 - x*y*z - y*z^2 + y*z);"
842  "4*( - x*y*z^2 + x*y*z);"
843  " - 2*x^2*y*z - 2*x*y^2*z + 2*x*y*z^2"
844  " + 2*x^2*z + 2*y^2*z - 2*x*z^2 - 2*y*z^2 + 3*x*y*z - x*z - y*z + 2*z^2 - z;"
845  "4*(x^2*y*z - x^2*z - x*y*z + x*z);"
846  " - 2*x^2*y*z + 2*x*y^2*z - 2*x*y*z^2 + 2*x^2*z + 2*x*z^2 + x*y*z - 3*x*z;"
847  "4*(x*y^2*z - y^2*z - x*y*z + y*z);"
848  "4*( - x*y^2*z + x*y*z);"
849  "2*x^2*y*z - 2*x*y^2*z - 2*x*y*z^2 + 2*y^2*z + 2*y*z^2 + x*y*z - 3*y*z;"
850  "4*( - x^2*y*z + x*y*z);"
851  "2*x^2*y*z + 2*x*y^2*z + 2*x*y*z^2 - 5*x*y*z;");
852 
853  for (int i = 0; i < 20; ++i)
854  trans[i] = read_base_poly(3, s);
855  }
856  fill_standard_vertices();
857  }
858  };
859 
860  static pgeometric_trans
861  Q2_incomplete_gt(gt_param_list& params,
862  std::vector<dal::pstatic_stored_object> &dependencies) {
863  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
864  << params.size() << " should be 1.");
865  GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
866  int n = int(::floor(params[0].num() + 0.01));
867  GMM_ASSERT1(n == 2 || n == 3, "Bad parameter, expected value 2 or 3");
868 
869  dependencies.push_back(Q2_incomplete_of_reference(dim_type(n)));
870  return std::make_shared<Q2_incomplete_trans_>(dim_type(n));
871  }
872 
873  pgeometric_trans Q2_incomplete_geotrans(dim_type nc) {
874  std::stringstream name;
875  name << "GT_Q2_INCOMPLETE(" << nc << ")";
876  return geometric_trans_descriptor(name.str());
877  }
878 
879  /* ******************************************************************** */
880  /* Pyramidal geometric transformation of order k=1 or 2. */
881  /* ******************************************************************** */
882 
883  struct pyramid_QK_trans_: public fraction_geometric_trans {
884  pyramid_QK_trans_(short_type k) {
885  cvr = pyramid_QK_of_reference(k);
886  size_type R = cvr->structure()->nb_points();
887  is_lin = false;
888  complexity_ = k;
889  trans.resize(R);
890 
891  if (k == 1) {
892  base_rational_fraction Q(read_base_poly(3, "x*y"), // Q = xy/(1-z)
893  read_base_poly(3, "1-z"));
894  trans[0] = (read_base_poly(3, "1-x-y-z") + Q)*0.25;
895  trans[1] = (read_base_poly(3, "1+x-y-z") - Q)*0.25;
896  trans[2] = (read_base_poly(3, "1-x+y-z") - Q)*0.25;
897  trans[3] = (read_base_poly(3, "1+x+y-z") + Q)*0.25;
898  trans[4] = read_base_poly(3, "z");
899  } else if (k == 2) {
900  base_poly xi0 = read_base_poly(3, "(1-z-x)*0.5");
901  base_poly xi1 = read_base_poly(3, "(1-z-y)*0.5");
902  base_poly xi2 = read_base_poly(3, "(1-z+x)*0.5");
903  base_poly xi3 = read_base_poly(3, "(1-z+y)*0.5");
904  base_poly x = read_base_poly(3, "x");
905  base_poly y = read_base_poly(3, "y");
906  base_poly z = read_base_poly(3, "z");
907  base_poly ones = read_base_poly(3, "1");
908  base_poly un_z = read_base_poly(3, "1-z");
909  base_rational_fraction Q(read_base_poly(3, "1"), un_z); // Q = 1/(1-z)
910  trans[ 0] = Q*Q*xi0*xi1*(x*y-z*un_z);
911  trans[ 1] = Q*Q*xi0*xi1*xi2*(xi1*2.-un_z)*4.;
912  trans[ 2] = Q*Q*xi1*xi2*(-x*y-z*un_z);
913  trans[ 3] = Q*Q*xi3*xi0*xi1*(xi0*2.-un_z)*4.;
914  trans[ 4] = Q*Q*xi0*xi1*xi2*xi3*16.;
915  trans[ 5] = Q*Q*xi1*xi2*xi3*(xi2*2.-un_z)*4.;
916  trans[ 6] = Q*Q*xi3*xi0*(-x*y-z*un_z);
917  trans[ 7] = Q*Q*xi2*xi3*xi0*(xi3*2.-un_z)*4.;
918  trans[ 8] = Q*Q*xi2*xi3*(x*y-z*un_z);
919  trans[ 9] = Q*z*xi0*xi1*4.;
920  trans[10] = Q*z*xi1*xi2*4.;
921  trans[11] = Q*z*xi3*xi0*4.;
922  trans[12] = Q*z*xi2*xi3*4.;
923  trans[13] = read_base_poly(3, "z*(2*z-1)");
924  }
925  fill_standard_vertices();
926  }
927  };
928 
929  static pgeometric_trans
930  pyramid_QK_gt(gt_param_list& params,
931  std::vector<dal::pstatic_stored_object> &deps) {
932  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
933  << params.size() << " should be 1.");
934  GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
935  int k = int(::floor(params[0].num() + 0.01));
936 
937  deps.push_back(pyramid_QK_of_reference(dim_type(k)));
938  return std::make_shared<pyramid_QK_trans_>(dim_type(k));
939  }
940 
941  pgeometric_trans pyramid_QK_geotrans(short_type k) {
942  static short_type k_ = -1;
943  static pgeometric_trans pgt = 0;
944  if (k != k_) {
945  std::stringstream name;
946  name << "GT_PYRAMID(" << k << ")";
947  pgt = geometric_trans_descriptor(name.str());
948  }
949  return pgt;
950  }
951 
952  /* ******************************************************************** */
953  /* Incomplete quadratic pyramidal geometric transformation. */
954  /* ******************************************************************** */
955 
956  struct pyramid_Q2_incomplete_trans_: public fraction_geometric_trans {
957  pyramid_Q2_incomplete_trans_() {
959  size_type R = cvr->structure()->nb_points();
960  is_lin = false;
961  complexity_ = 2;
962  trans.resize(R);
963 
964  base_poly xy = read_base_poly(3, "x*y");
965  base_poly z = read_base_poly(3, "z");
966 
967  base_poly p00m = read_base_poly(3, "1-z");
968 
969  base_poly pmmm = read_base_poly(3, "1-x-y-z");
970  base_poly ppmm = read_base_poly(3, "1+x-y-z");
971  base_poly pmpm = read_base_poly(3, "1-x+y-z");
972  base_poly pppm = read_base_poly(3, "1+x+y-z");
973 
974  base_poly mmm0_ = read_base_poly(3, "-1-x-y")*0.25;
975  base_poly mpm0_ = read_base_poly(3, "-1+x-y")*0.25;
976  base_poly mmp0_ = read_base_poly(3, "-1-x+y")*0.25;
977  base_poly mpp0_ = read_base_poly(3, "-1+x+y")*0.25;
978 
979  base_poly pp0m = read_base_poly(3, "1+x-z");
980  base_poly pm0m = read_base_poly(3, "1-x-z");
981  base_poly p0mm = read_base_poly(3, "1-y-z");
982  base_poly p0pm = read_base_poly(3, "1+y-z");
983 
984  trans[ 0] = mmm0_ * pmmm + base_rational_fraction(mmm0_ * xy, p00m); // N5 (Bedrosian, 1992)
985  trans[ 1] = base_rational_fraction(pp0m * pm0m * p0mm * 0.5, p00m); // N6
986  trans[ 2] = mpm0_ * ppmm - base_rational_fraction(mpm0_ * xy, p00m); // N7
987  trans[ 3] = base_rational_fraction(p0pm * p0mm * pm0m * 0.5, p00m); // N4
988  trans[ 4] = base_rational_fraction(p0pm * p0mm * pp0m * 0.5, p00m); // N8
989  trans[ 5] = mmp0_ * pmpm - base_rational_fraction(mmp0_ * xy, p00m); // N3
990  trans[ 6] = base_rational_fraction(pp0m * pm0m * p0pm * 0.5, p00m); // N2
991  trans[ 7] = mpp0_ * pppm + base_rational_fraction(mpp0_ * xy, p00m); // N1
992  trans[ 8] = base_rational_fraction(pm0m * p0mm * z, p00m); // N11
993  trans[ 9] = base_rational_fraction(pp0m * p0mm * z, p00m); // N12
994  trans[10] = base_rational_fraction(pm0m * p0pm * z, p00m); // N10
995  trans[11] = base_rational_fraction(pp0m * p0pm * z, p00m); // N9
996  trans[12] = read_base_poly(3, "2*z^2-z"); // N13
997 
998  fill_standard_vertices();
999  }
1000  };
1001 
1002  static pgeometric_trans
1003  pyramid_Q2_incomplete_gt(gt_param_list& params,
1004  std::vector<dal::pstatic_stored_object> &deps) {
1005  GMM_ASSERT1(params.size() == 0, "Bad number of parameters : "
1006  << params.size() << " should be 0.");
1007 
1008  deps.push_back(pyramid_Q2_incomplete_of_reference());
1009  return std::make_shared<pyramid_Q2_incomplete_trans_>();
1010  }
1011 
1012  pgeometric_trans pyramid_Q2_incomplete_geotrans() {
1013  static pgeometric_trans pgt = 0;
1014  if (!pgt)
1015  pgt = geometric_trans_descriptor("GT_PYRAMID_Q2_INCOMPLETE");
1016  return pgt;
1017  }
1018 
1019  /* ******************************************************************** */
1020  /* Incomplete quadratic prism geometric transformation. */
1021  /* ******************************************************************** */
1022 
1023  struct prism_incomplete_P2_trans_: public poly_geometric_trans {
1024  prism_incomplete_P2_trans_() {
1026  size_type R = cvr->structure()->nb_points();
1027  is_lin = false;
1028  complexity_ = 2;
1029  trans.resize(R);
1030 
1031  std::stringstream s
1032  ( "-2*y*z^2-2*x*z^2+2*z^2-2*y^2*z-4*x*y*z+5*y*z-2*x^2*z+5*x*z"
1033  "-3*z+2*y^2+4*x*y-3*y+2*x^2-3*x+1;"
1034  "4*(x*y*z+x^2*z-x*z-x*y-x^2+x);"
1035  "2*x*z^2-2*x^2*z-x*z+2*x^2-x;"
1036  "4*(y^2*z+x*y*z-y*z-y^2-x*y+y);"
1037  "4*(x*y-x*y*z);"
1038  "2*y*z^2-2*y^2*z-y*z+2*y^2-y;"
1039  "4*(y*z^2+x*z^2-z^2-y*z-x*z+z);"
1040  "4*(x*z-x*z^2);"
1041  "4*(y*z-y*z^2);"
1042  "-2*y*z^2-2*x*z^2+2*z^2+2*y^2*z+4*x*y*z-y*z+2*x^2*z-x*z-z;"
1043  "4*(-x*y*z-x^2*z+x*z);"
1044  "2*x*z^2+2*x^2*z-3*x*z;"
1045  "4*(-y^2*z-x*y*z+y*z);"
1046  "4*x*y*z;"
1047  "2*y*z^2+2*y^2*z-3*y*z;");
1048 
1049  for (int i = 0; i < 15; ++i)
1050  trans[i] = read_base_poly(3, s);
1051 
1052  fill_standard_vertices();
1053  }
1054  };
1055 
1056  static pgeometric_trans
1057  prism_incomplete_P2_gt(gt_param_list& params,
1058  std::vector<dal::pstatic_stored_object> &deps) {
1059  GMM_ASSERT1(params.size() == 0, "Bad number of parameters : "
1060  << params.size() << " should be 0.");
1061 
1062  deps.push_back(prism_incomplete_P2_of_reference());
1063  return std::make_shared<prism_incomplete_P2_trans_>();
1064  }
1065 
1066  pgeometric_trans prism_incomplete_P2_geotrans() {
1067  static pgeometric_trans pgt = 0;
1068  if (!pgt)
1069  pgt = geometric_trans_descriptor("GT_PRISM_INCOMPLETE_P2");
1070  return pgt;
1071  }
1072 
1073  /* ******************************************************************** */
1074  /* Misc function. */
1075  /* ******************************************************************** */
1076 
1077  /* norm of returned vector is the ratio between the face surface on
1078  the reference element and the face surface on the real element.
1079  IT IS NOT UNITARY
1080  */
1081  base_small_vector
1083  GMM_ASSERT1(c.G().ncols() == c.pgt()->nb_points(), "dimensions mismatch");
1084  base_small_vector un(c.N());
1085  gmm::mult(c.B(), c.pgt()->normals()[face], un);
1086  return un;
1087  }
1088 
1089  /*
1090  return the local basis (i.e. the normal in the first column, and the
1091  tangent vectors in the other columns
1092  */
1093  base_matrix
1095  size_type face) {
1096  GMM_ASSERT1(c.G().ncols() == c.pgt()->nb_points(), "dimensions mismatch");
1097  base_small_vector up = c.pgt()->normals()[face];
1098  size_type P = c.pgt()->structure()->dim();
1099 
1100  base_matrix baseP(P, P);
1101  gmm::copy(gmm::identity_matrix(), baseP);
1102  size_type i0 = 0;
1103  for (size_type i=1; i < P; ++i)
1104  if (gmm::abs(up[i])>gmm::abs(up[i0])) i0=i;
1105  if (i0) gmm::copy(gmm::mat_col(baseP, 0), gmm::mat_col(baseP, i0));
1106  gmm::copy(up, gmm::mat_col(baseP, 0));
1107 
1108  base_matrix baseN(c.N(), P);
1109  gmm::mult(c.B(), baseP, baseN);
1110 
1111  /* Modified Gram-Schmidt */
1112  for (size_type k=0; k < P; ++k) {
1113  for (size_type l=0; l < k; ++l) {
1114  gmm::add(gmm::scaled(gmm::mat_col(baseN,l),
1115  -gmm::vect_sp(gmm::mat_col(baseN,l),
1116  gmm::mat_col(baseN,k))),
1117  gmm::mat_col(baseN,k));
1118  }
1119  gmm::scale(gmm::mat_col(baseN,k),
1120  1./gmm::vect_norm2(gmm::mat_col(baseN,k)));
1121  }
1122  /* TODO: for cases where P < N, complete the basis */
1123 
1124  /* Ensure that the baseN is direct */
1125  if (c.N() == P && c.N()>1 && gmm::lu_det(baseN) < 0) {
1126  gmm::scale(gmm::mat_col(baseN,1),-1.);
1127  }
1128  return baseN;
1129  }
1130 
1131 
1132 
1133 
1134  /* ******************************************************************** */
1135  /* Naming system */
1136  /* ******************************************************************** */
1137 
1138  struct geometric_trans_naming_system
1139  : public dal::naming_system<geometric_trans> {
1140  geometric_trans_naming_system() :
1141  dal::naming_system<geometric_trans>("GT") {
1142  add_suffix("PK", PK_gt);
1143  add_suffix("QK", QK_gt);
1144  add_suffix("PRISM_PK", prism_pk_gt);
1145  add_suffix("PRISM", prism_pk_gt);
1146  add_suffix("PRODUCT", product_gt);
1147  add_suffix("LINEAR_PRODUCT", linear_product_gt);
1148  add_suffix("LINEAR_QK", linear_qk);
1149  add_suffix("Q2_INCOMPLETE", Q2_incomplete_gt);
1150  add_suffix("PYRAMID_QK", pyramid_QK_gt);
1151  add_suffix("PYRAMID", pyramid_QK_gt);
1152  add_suffix("PYRAMID_Q2_INCOMPLETE", pyramid_Q2_incomplete_gt);
1153  add_suffix("PRISM_INCOMPLETE_P2", prism_incomplete_P2_gt);
1154  }
1155  };
1156 
1157  void add_geometric_trans_name
1158  (std::string name, dal::naming_system<geometric_trans>::pfunction f) {
1160  f);
1161  }
1162 
1164  size_type i = 0;
1166  auto ptrans = name_system.method(name, i);
1167  auto &trans = const_cast<bgeot::geometric_trans&>(*ptrans);
1168  auto short_name = name_system.shorter_name_of_method(ptrans);
1169  trans.set_name(short_name);
1170  return ptrans;
1171  }
1172 
1175  const torus_geom_trans *pgt_torus = dynamic_cast<const torus_geom_trans *>(p.get());
1176  if (pgt_torus) return instance.shorter_name_of_method(pgt_torus->get_original_transformation());
1177  return instance.shorter_name_of_method(p);
1178  }
1179 
1180  /* Fonctions pour la ref. directe. */
1181 
1182  pgeometric_trans simplex_geotrans(size_type n, short_type k) {
1183  static pgeometric_trans pgt = 0;
1184  static size_type d = size_type(-2);
1185  static short_type r = short_type(-2);
1186  if (d != n || r != k) {
1187  std::stringstream name;
1188  name << "GT_PK(" << n << "," << k << ")";
1189  pgt = geometric_trans_descriptor(name.str());
1190  d = n; r = k;
1191  }
1192  return pgt;
1193  }
1194 
1195  pgeometric_trans parallelepiped_geotrans(size_type n, short_type k) {
1196  static pgeometric_trans pgt = 0;
1197  static size_type d = size_type(-2);
1198  static short_type r = short_type(-2);
1199  if (d != n || r != k) {
1200  std::stringstream name;
1201  name << "GT_QK(" << n << "," << k << ")";
1202  pgt = geometric_trans_descriptor(name.str());
1203  d = n; r = k;
1204  }
1205  return pgt;
1206  }
1207 
1208  static std::string name_of_linear_qk_trans(size_type dim) {
1209  switch (dim) {
1210  case 1: return "GT_PK(1,1)";
1211  default: return std::string("GT_LINEAR_PRODUCT(")
1212  + name_of_linear_qk_trans(dim-1)
1213  + std::string(",GT_PK(1,1))");
1214  }
1215  }
1216 
1217  pgeometric_trans parallelepiped_linear_geotrans(size_type n) {
1218  static pgeometric_trans pgt = 0;
1219  static size_type d = size_type(-2);
1220  if (d != n) {
1221  std::stringstream name(name_of_linear_qk_trans(n));
1222  pgt = geometric_trans_descriptor(name.str());
1223  d = n;
1224  }
1225  return pgt;
1226  }
1227 
1228  pgeometric_trans prism_linear_geotrans(size_type n) {
1229  static pgeometric_trans pgt = 0;
1230  static size_type d = size_type(-2);
1231  if (d != n) {
1232  std::stringstream name;
1233  name << "GT_LINEAR_PRODUCT(GT_PK(" << (n-1) << ", 1), GT_PK(1,1))";
1234  pgt = geometric_trans_descriptor(name.str());
1235  d = n;
1236  }
1237  return pgt;
1238  }
1239 
1240  pgeometric_trans linear_product_geotrans(pgeometric_trans pg1,
1241  pgeometric_trans pg2) {
1242  std::stringstream name;
1243  name << "GT_LINEAR_PRODUCT(" << name_of_geometric_trans(pg1) << ","
1244  << name_of_geometric_trans(pg2) << ")";
1245  return geometric_trans_descriptor(name.str());
1246  }
1247 
1248  pgeometric_trans prism_geotrans(size_type n, short_type k) {
1249  static pgeometric_trans pgt = 0;
1250  static size_type d = size_type(-2);
1251  static short_type r = short_type(-2);
1252  if (d != n || r != k) {
1253  std::stringstream name;
1254  name << "GT_PRISM(" << n << "," << k << ")";
1255  pgt = geometric_trans_descriptor(name.str());
1256  d = n; r = k;
1257  }
1258  return pgt;
1259  }
1260 
1261  pgeometric_trans product_geotrans(pgeometric_trans pg1,
1262  pgeometric_trans pg2) {
1263  static pgeometric_trans pgt = 0;
1264  static pgeometric_trans pg1_ = 0;
1265  static pgeometric_trans pg2_ = 0;
1266  if (pg1 != pg1_ || pg2 != pg2_) {
1267  std::stringstream name;
1268  name << "GT_PRODUCT(" << name_of_geometric_trans(pg1) << ","
1269  << name_of_geometric_trans(pg2) << ")";
1270  pgt = geometric_trans_descriptor(name.str());
1271  pg1_ = pg1; pg2_ = pg2;
1272  }
1273  return pgt;
1274  }
1275 
1276 
1277  // To be completed
1278  pgeometric_trans default_trans_of_cvs(pconvex_structure cvs) {
1279 
1280  dim_type n = cvs->dim();
1281  short_type nbf = cvs->nb_faces();
1282  short_type nbpt = cvs->nb_points();
1283 
1284  // Basic cases
1285  if (cvs == simplex_structure(n)) return simplex_geotrans(n, 1);
1286  if (cvs == parallelepiped_structure(n))
1287  return parallelepiped_geotrans(n, 1);
1288  if (cvs == prism_P1_structure(n)) return prism_geotrans(n, 1);
1289 
1290  // more elaborated ones
1291  switch (n) {
1292  case 1 : return simplex_geotrans(1, short_type(nbpt-1));
1293  case 2 :
1294  if (nbf == 3) {
1295  short_type k = short_type(round((sqrt(1.+8.*nbpt) - 3. ) / 2.));
1296  if (cvs == simplex_structure(2,k)) return simplex_geotrans(2, k);
1297  } else if (nbf == 4) {
1298  short_type k = short_type(round(sqrt(1.*nbpt)) - 1.);
1299  if (cvs == parallelepiped_structure(2, k))
1300  return parallelepiped_geotrans(2, k);
1301  }
1302  break;
1303  case 3 :
1304  if (nbf == 4) {
1305  if (cvs == simplex_structure(3, 2)) return simplex_geotrans(3, 2);
1306  if (cvs == simplex_structure(3, 3)) return simplex_geotrans(3, 3);
1307  if (cvs == simplex_structure(3, 4)) return simplex_geotrans(3, 4);
1308  if (cvs == simplex_structure(3, 5)) return simplex_geotrans(3, 5);
1309  if (cvs == simplex_structure(3, 6)) return simplex_geotrans(3, 6);
1310  } else if (nbf == 6) {
1311  short_type k = short_type(round(pow(1.*nbpt, 1./3.)) - 1.);
1312  if (cvs == parallelepiped_structure(3, k))
1313  return parallelepiped_geotrans(3, k);
1314  } else if (nbf == 5) {
1315  if (cvs == pyramid_QK_structure(1)) return pyramid_QK_geotrans(1);
1316  if (cvs == pyramid_QK_structure(2)) return pyramid_QK_geotrans(2);
1317  if (cvs == pyramid_QK_structure(3)) return pyramid_QK_geotrans(3);
1318  if (cvs == pyramid_QK_structure(4)) return pyramid_QK_geotrans(4);
1319  if (cvs == pyramid_QK_structure(5)) return pyramid_QK_geotrans(5);
1320  if (cvs == pyramid_QK_structure(6)) return pyramid_QK_geotrans(6);
1321  if (cvs == pyramid_Q2_incomplete_structure())
1322  return pyramid_Q2_incomplete_geotrans();
1323  }
1324  break;
1325  }
1326  GMM_ASSERT1(false, "Unrecognized structure");
1327  }
1328 
1329 
1330 
1331  /* ********************************************************************* */
1332  /* Precomputation on geometric transformations. */
1333  /* ********************************************************************* */
1334 
1335  DAL_DOUBLE_KEY(pre_geot_key_, pgeometric_trans, pstored_point_tab);
1336 
1337  geotrans_precomp_::geotrans_precomp_(pgeometric_trans pg,
1338  pstored_point_tab ps)
1339  : pgt(pg), pspt(ps)
1340  { DAL_STORED_OBJECT_DEBUG_CREATED(this, "Geotrans precomp"); }
1341 
1342  void geotrans_precomp_::init_val() const {
1343  c.clear();
1344  c.resize(pspt->size(), base_vector(pgt->nb_points()));
1345  for (size_type j = 0; j < pspt->size(); ++j)
1346  pgt->poly_vector_val((*pspt)[j], c[j]);
1347  }
1348 
1349  void geotrans_precomp_::init_grad() const {
1350  dim_type N = pgt->dim();
1351  pc.clear();
1352  pc.resize(pspt->size(), base_matrix(pgt->nb_points() , N));
1353  for (size_type j = 0; j < pspt->size(); ++j)
1354  pgt->poly_vector_grad((*pspt)[j], pc[j]);
1355  }
1356 
1357  void geotrans_precomp_::init_hess() const {
1358  base_poly P, Q;
1359  dim_type N = pgt->structure()->dim();
1360  hpc.clear();
1361  hpc.resize(pspt->size(), base_matrix(pgt->nb_points(), gmm::sqr(N)));
1362  for (size_type j = 0; j < pspt->size(); ++j)
1363  pgt->poly_vector_hess((*pspt)[j], hpc[j]);
1364  }
1365 
1366  base_node geotrans_precomp_::transform(size_type i,
1367  const base_matrix &G) const {
1368  if (c.empty()) init_val();
1369  size_type N = G.nrows(), k = pgt->nb_points();
1370  base_node P(N);
1371  base_matrix::const_iterator git = G.begin();
1372  for (size_type l = 0; l < k; ++l) {
1373  scalar_type a = c[i][l];
1374  base_node::iterator pit = P.begin(), pite = P.end();
1375  for (; pit != pite; ++git, ++pit) *pit += a * (*git);
1376  }
1377  return P;
1378  }
1379 
1380  pgeotrans_precomp geotrans_precomp(pgeometric_trans pg,
1381  pstored_point_tab pspt,
1382  dal::pstatic_stored_object dep) {
1383  dal::pstatic_stored_object_key pk= std::make_shared<pre_geot_key_>(pg,pspt);
1384  dal::pstatic_stored_object o = dal::search_stored_object(pk);
1385  if (o) return std::dynamic_pointer_cast<const geotrans_precomp_>(o);
1386  pgeotrans_precomp p = std::make_shared<geotrans_precomp_>(pg, pspt);
1387  dal::add_stored_object(pk, p, pg, pspt, dal::AUTODELETE_STATIC_OBJECT);
1388  if (dep) dal::add_dependency(p, dep);
1389  return p;
1390  }
1391 
1392  void delete_geotrans_precomp(pgeotrans_precomp pgp)
1393  { dal::del_stored_object(pgp, true); }
1394 
1395 } /* end of namespace bgeot. */
1396 
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:165
bgeot::geotrans_interpolation_context::cv_center
const base_node & cv_center() const
coordinates of the center of the real convex.
Definition: bgeot_geometric_trans.cc:308
bgeot::geometric_trans
Description of a geometric transformation between a reference element and a real element.
Definition: bgeot_geometric_trans.h:105
bgeot::pyramid_Q2_incomplete_of_reference
pconvex_ref pyramid_Q2_incomplete_of_reference()
incomplete quadratic pyramidal element of reference (13-node)
Definition: bgeot_convex_ref.cc:576
bgeot::geotrans_interpolation_context::G
const base_matrix & G() const
matrix whose columns are the vertices of the convex
Definition: bgeot_geometric_trans.h:451
bgeot::compute_normal
base_small_vector compute_normal(const geotrans_interpolation_context &c, size_type face)
norm of returned vector is the ratio between the face surface on the real element and the face surfac...
Definition: bgeot_geometric_trans.cc:1082
bgeot::name_of_geometric_trans
std::string name_of_geometric_trans(pgeometric_trans p)
Get the string name of a geometric transformation.
Definition: bgeot_geometric_trans.cc:1173
bgeot_torus.h
Provides mesh of torus.
bgeot::geotrans_interpolation_context::xreal
const base_node & xreal() const
coordinates of the current point, in the real convex.
Definition: bgeot_geometric_trans.cc:299
bgeot::prism_P1_structure
pconvex_structure prism_P1_structure(dim_type nc)
Give a pointer on the structures of a prism of dimension d.
Definition: bgeot_convex_structure.h:206
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
dal_singleton.h
A simple singleton implementation.
bgeot::geotrans_interpolation_context::K
const base_matrix & K() const
See getfem kernel doc for these matrices.
Definition: bgeot_geometric_trans.cc:358
bgeot::pyramid_QK_of_reference
pconvex_ref pyramid_QK_of_reference(dim_type k)
pyramidal element of reference of degree k (k = 1 or 2 only)
Definition: bgeot_convex_ref.cc:522
dal
Dynamic Array Library.
Definition: dal_backtrace.cc:29
bgeot::geometric_trans::poly_vector_val
virtual void poly_vector_val(const base_node &pt, base_vector &val) const =0
Gives the value of the functions vector at a certain point.
dal::naming_system
Associate a name to a method descriptor and store method descriptors.
Definition: dal_naming_system.h:56
bgeot::geotrans_interpolation_context
the geotrans_interpolation_context structure is passed as the argument of geometric transformation in...
Definition: bgeot_geometric_trans.h:411
dal::del_stored_object
void del_stored_object(const pstatic_stored_object &o, bool ignore_unstored)
Delete an object and the object which depend on it.
Definition: dal_static_stored_objects.cc:369
bgeot::geotrans_interpolation_context::G_
const base_matrix * G_
transformed point
Definition: bgeot_geometric_trans.h:415
bgeot::geometric_trans_descriptor
pgeometric_trans geometric_trans_descriptor(std::string name)
Get the geometric transformation from its string name.
Definition: bgeot_geometric_trans.cc:1163
bgeot::geotrans_interpolation_context::set_xref
void set_xref(const base_node &P)
change the current point (coordinates given in the reference convex)
Definition: bgeot_geometric_trans.cc:458
bgeot::geotrans_interpolation_context::xref
const base_node & xref() const
coordinates of the current point, in the reference convex.
Definition: bgeot_geometric_trans.cc:291
bgeot::geometric_trans::dim
dim_type dim() const
Dimension of the reference element.
Definition: bgeot_geometric_trans.h:119
bgeot::read_base_poly
base_poly read_base_poly(short_type n, std::istream &f)
read a base_poly on the stream ist.
Definition: bgeot_poly.cc:211
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::simplex_structure
pconvex_structure simplex_structure(dim_type nc)
Give a pointer on the structures of a simplex of dimension d.
Definition: bgeot_convex_structure.cc:148
dal_tree_sorted.h
a balanced tree stored in a dal::dynamic_array
dal::search_stored_object
pstatic_stored_object search_stored_object(pstatic_stored_object_key k)
Gives a pointer to an object from a key pointer.
Definition: dal_static_stored_objects.cc:177
bgeot::convex_ref_product
pconvex_ref convex_ref_product(pconvex_ref a, pconvex_ref b)
tensorial product of two convex ref.
Definition: bgeot_convex_ref.cc:729
bgeot::prism_incomplete_P2_of_reference
pconvex_ref prism_incomplete_P2_of_reference()
incomplete quadratic prism element of reference (15-node)
Definition: bgeot_convex_ref.cc:635
bgeot::geotrans_interpolation_context::K_
base_matrix K_
real center of convex (average of columns of G)
Definition: bgeot_geometric_trans.h:417
bgeot::geotrans_interpolation_context::cv_center_
base_node cv_center_
pointer to the matrix of real nodes of the convex
Definition: bgeot_geometric_trans.h:416
bgeot::geotrans_interpolation_context::xreal_
base_node xreal_
reference point
Definition: bgeot_geometric_trans.h:414
bgeot::small_vector< scalar_type >
bgeot_geometric_trans.h
Geometric transformations on convexes.
bgeot::geotrans_interpolation_context::PC
base_matrix PC
Jacobian.
Definition: bgeot_geometric_trans.h:423
bgeot::pyramid_Q2_incomplete_structure
pconvex_structure pyramid_Q2_incomplete_structure()
Give a pointer on the 3D quadratic incomplete pyramid structure.
Definition: bgeot_convex_structure.cc:605
bgeot::geotrans_interpolation_context::J_
scalar_type J_
index of current point in the pgp
Definition: bgeot_geometric_trans.h:422
bgeot::geotrans_interpolation_context::pgt_
pgeometric_trans pgt_
see documentation (getfem kernel doc) for more details
Definition: bgeot_geometric_trans.h:418
dal::add_stored_object
void add_stored_object(pstatic_stored_object_key k, pstatic_stored_object o, permanence perm)
Add an object with two optional dependencies.
Definition: dal_static_stored_objects.cc:284
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:27
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:54
bgeot::torus_geom_trans
An adaptor that adapts a two dimensional geometric_trans to include radial dimension.
Definition: bgeot_torus.h:48
bgeot::geometric_trans::compute_K_matrix
virtual void compute_K_matrix(const base_matrix &G, const base_matrix &pc, base_matrix &K) const
compute K matrix from multiplication of G with gradient
Definition: bgeot_geometric_trans.cc:272
bgeot::parallelepiped_structure
pconvex_structure parallelepiped_structure(dim_type nc, dim_type k)
Give a pointer on the structures of a parallelepiped of dimension d.
Definition: bgeot_convex_structure.cc:397
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
bgeot::simplex_of_reference
pconvex_ref simplex_of_reference(dim_type nc, short_type K)
returns a simplex of reference of dimension nc and degree k
Definition: bgeot_convex_ref.cc:340
bgeot::pyramid_QK_structure
pconvex_structure pyramid_QK_structure(dim_type k)
Give a pointer on the 3D pyramid structure for a degree k = 1 or 2.
Definition: bgeot_convex_structure.cc:508
bgeot::geometric_trans::nb_points
size_type nb_points() const
Number of geometric nodes.
Definition: bgeot_geometric_trans.h:123
bgeot_poly_composite.h
Handle composite polynomials.
bgeot::geotrans_interpolation_context::ii_
size_type ii_
if pgp != 0, it is the same as pgp's one
Definition: bgeot_geometric_trans.h:421
bgeot::geometric_trans::transform
base_node transform(const base_node &pt, const CONT &PTAB) const
Apply the geometric transformation to point pt, PTAB contains the points of the real convex.
Definition: bgeot_geometric_trans.h:175
bgeot::Q2_incomplete_of_reference
pconvex_ref Q2_incomplete_of_reference(dim_type nc)
incomplete Q2 quadrilateral/hexahedral of reference of dimension d = 2 or 3
Definition: bgeot_convex_ref.cc:426
dal::add_dependency
void add_dependency(pstatic_stored_object o1, pstatic_stored_object o2)
Add a dependency, object o1 will depend on object o2.
Definition: dal_static_stored_objects.cc:230
bgeot::compute_local_basis
base_matrix compute_local_basis(const geotrans_interpolation_context &c, size_type face)
return the local basis (i.e.
Definition: bgeot_geometric_trans.cc:1094