GetFEM  5.4.2
getfem_mat_elem.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 <deque>
24 #include "getfem/dal_singleton.h"
25 #include "getfem/getfem_fem.h"
26 #include "getfem/getfem_mat_elem.h"
27 #include "getfem/getfem_omp.h"
28 
29 namespace getfem {
30  /* ********************************************************************* */
31  /* Elementary matrices computation. */
32  /* ********************************************************************* */
33 
34  struct emelem_comp_key_ : virtual public dal::static_stored_object_key {
35  pmat_elem_type pmt;
36  pintegration_method ppi;
37  bgeot::pgeometric_trans pgt;
38  /* prefer_comp_on_real_element: compute elementary matrices on the real
39  element if possible (i.e. if no exact integration is used); this allow
40  using inline reduction during the integration */
41  bool prefer_comp_on_real_element;
42  bool compare(const static_stored_object_key &oo) const override{
43  auto &o = dynamic_cast<const emelem_comp_key_ &>(oo);
44  if (pmt < o.pmt) return true;
45  if (o.pmt < pmt) return false;
46  if (ppi < o.ppi) return true;
47  if (o.ppi < ppi) return false;
48  if (pgt < o.pgt) return true;
49  if (o.pgt < pgt) return false;
50  return prefer_comp_on_real_element < o.prefer_comp_on_real_element;
51  }
52  bool equal(const static_stored_object_key &oo) const override{
53  auto &o = dynamic_cast<const emelem_comp_key_ &>(oo);
54 
55  if (pmt == o.pmt && ppi == o.ppi && pgt == o.pgt) return true;
56 
57  auto pmat_key = dal::key_of_stored_object(pmt);
58  auto poo_mat_key = dal::key_of_stored_object(o.pmt);
59  if (*pmat_key != *poo_mat_key) return false;
60 
61  auto pint_key = dal::key_of_stored_object(ppi);
62  auto poo_int_key = dal::key_of_stored_object(o.ppi);
63  if (*pint_key != *poo_int_key) return false;
64 
65  auto pgt_key = dal::key_of_stored_object(pgt);
66  auto poo_gt_key = dal::key_of_stored_object(o.pgt);
67  if (*pgt_key != *poo_gt_key) return false;
68 
69  return true;
70  }
71  emelem_comp_key_(pmat_elem_type pm, pintegration_method pi,
72  bgeot::pgeometric_trans pg, bool on_relt)
73  { pmt = pm; ppi = pi; pgt = pg; prefer_comp_on_real_element = on_relt; }
74  emelem_comp_key_(void) { }
75  };
76 
77  struct emelem_comp_structure_ : public mat_elem_computation {
78  bgeot::pgeotrans_precomp pgp;
79  ppoly_integration ppi;
80  papprox_integration pai;
81  bool is_ppi;
82  mutable std::vector<base_tensor> mref;
83  mutable std::vector<pfem_precomp> pfp;
84  mutable std::vector<base_tensor> elmt_stored;
85  short_type nbf, dim;
86  std::deque<short_type> grad_reduction, hess_reduction, trans_reduction;
87  std::deque<short_type> K_reduction;
88  std::deque<pfem> trans_reduction_pfi;
89  mutable base_vector un, up;
90  mutable bool faces_computed;
91  mutable bool volume_computed;
92  bool is_linear;
93  bool computed_on_real_element;
94  size_type memsize() const {
95  size_type sz = sizeof(emelem_comp_structure_) +
96  mref.capacity()*sizeof(base_tensor) +
97  grad_reduction.size()*sizeof(short_type) +
98  K_reduction.size()*sizeof(short_type) +
99  hess_reduction.size()*sizeof(short_type) +
100  trans_reduction.size()*sizeof(short_type) +
101  trans_reduction_pfi.size()*sizeof(pfem);
102 
103  for (size_type i=0; i < mref.size(); ++i) sz += mref[i].memsize();
104  return sz;
105  }
106 
107  emelem_comp_structure_(pmat_elem_type pm, pintegration_method pi,
108  bgeot::pgeometric_trans pg,
109  bool prefer_comp_on_real_element) {
110 
111  pgt = pg;
112  pgp = bgeot::geotrans_precomp(pg, pi->pintegration_points(), pi);
113  pme = pm;
114  switch (pi->type()) {
115  case IM_EXACT:
116  ppi = pi->exact_method(); pai = 0; is_ppi = true; break;
117  case IM_APPROX:
118  ppi = 0; pai = pi->approx_method(); is_ppi = false; break;
119  case IM_NONE:
120  GMM_ASSERT1(false, "Attempt to use IM_NONE integration method "
121  "in assembly!\n");
122  }
123 
124  faces_computed = volume_computed = false;
125  is_linear = pgt->is_linear();
126  computed_on_real_element = !is_linear || (prefer_comp_on_real_element && !is_ppi);
127  // computed_on_real_element = true;
128  nbf = pgt->structure()->nb_faces();
129  dim = pgt->structure()->dim();
130  mat_elem_type::const_iterator it = pme->begin(), ite = pme->end();
131  // size_type d = pgt->dim();
132 
133  for (short_type k = 0; it != ite; ++it, ++k) {
134  if ((*it).pfi) {
135  if ((*it).pfi->is_on_real_element()) computed_on_real_element = true;
136  GMM_ASSERT1(!is_ppi || (((*it).pfi->is_polynomial()) && is_linear
137  && !computed_on_real_element),
138  "Exact integration not allowed in this context");
139 
140  if ((*it).t != GETFEM_NONLINEAR_ && !((*it).pfi->is_equivalent())) {
141  // TODO : le numero d'indice � reduire peut changer ...
142  trans_reduction.push_back(k);
143  trans_reduction_pfi.push_back((*it).pfi);
144  }
145  }
146  switch ((*it).t) {
147  case GETFEM_BASE_ :
148  if ((*it).pfi->target_dim() > 1) {
149  ++k;
150  switch((*it).pfi->vectorial_type()) {
151  case virtual_fem::VECTORIAL_PRIMAL_TYPE:
152  K_reduction.push_back(k); break;
153  case virtual_fem::VECTORIAL_DUAL_TYPE:
154  grad_reduction.push_back(k); break;
155  default: break;
156  }
157  }
158  break;
159  case GETFEM_UNIT_NORMAL_ :
160  computed_on_real_element = true; break;
161  case GETFEM_GRAD_GEOTRANS_ :
162  case GETFEM_GRAD_GEOTRANS_INV_ :
163  ++k; computed_on_real_element = true; break;
164  case GETFEM_GRAD_ : {
165  ++k;
166  switch((*it).pfi->vectorial_type()) {
167  case virtual_fem::VECTORIAL_PRIMAL_TYPE:
168  K_reduction.push_back(k); break;
169  case virtual_fem::VECTORIAL_DUAL_TYPE:
170  grad_reduction.push_back(k); break;
171  default: break;
172  }
173  if ((*it).pfi->target_dim() > 1) ++k;
174  if (!((*it).pfi->is_on_real_element()))
175  grad_reduction.push_back(k);
176  } break;
177  case GETFEM_HESSIAN_ : {
178  ++k;
179  switch((*it).pfi->vectorial_type()) {
180  case virtual_fem::VECTORIAL_PRIMAL_TYPE:
181  K_reduction.push_back(k); break;
182  case virtual_fem::VECTORIAL_DUAL_TYPE:
183  grad_reduction.push_back(k); break;
184  default: break;
185  }
186 
187  if ((*it).pfi->target_dim() > 1) ++k;
188  if (!((*it).pfi->is_on_real_element()))
189  hess_reduction.push_back(k);
190  } break;
191  case GETFEM_NONLINEAR_ : {
192  if ((*it).nl_part == 0) {
193  k = short_type(k+(*it).nlt->sizes(size_type(-1)).size()-1);
194  GMM_ASSERT1(!is_ppi, "For nonlinear terms you have "
195  "to use approximated integration");
196  computed_on_real_element = true;
197  }
198  } break;
199  }
200  }
201 
202  if (!is_ppi) {
203  pfp.resize(pme->size());
204  it = pme->begin(), ite = pme->end();
205  for (size_type k = 0; it != ite; ++it, ++k)
206  if ((*it).pfi)
207  pfp[k] = fem_precomp((*it).pfi, pai->pintegration_points(), pi);
208  else pfp[k] = 0;
209  elmt_stored.resize(pme->size());
210  }
211  if (!computed_on_real_element) mref.resize(nbf + 1);
212  }
213 
214  void add_elem(base_tensor &t, fem_interpolation_context& ctx,
215  scalar_type J, bool first, bool trans,
216  mat_elem_integration_callback *icb,
217  bgeot::multi_index sizes) const {
218  mat_elem_type::const_iterator it = pme->begin(), ite = pme->end();
219  bgeot::multi_index aux_ind;
220 
221  for (size_type k = 0; it != ite; ++it, ++k) {
222  if ((*it).t == GETFEM_NONLINEAR_)
223  (*it).nlt->term_num() = size_type(-1);
224  }
225  it = pme->begin();
226 
227  // incrementing "mit" should match increments of "j" in mat_elem_type::sizes
228  bgeot::multi_index::iterator mit = sizes.begin();
229  for (size_type k = 0; it != ite; ++it, ++k, ++mit) {
230  if (pfp[k]) ctx.set_pfp(pfp[k]);
231 
232  switch ((*it).t) {
233  case GETFEM_BASE_ :
234  if ((*it).pfi && (*it).pfi->target_dim() > 1) ++mit;
235  if (trans)
236  (*it).pfi->real_base_value(ctx, elmt_stored[k], icb != 0);
237  else
238  elmt_stored[k] = pfp[k]->val(ctx.ii());
239  break;
240  case GETFEM_GRAD_ :
241  ++mit;
242  if ((*it).pfi && (*it).pfi->target_dim() > 1) ++mit;
243  if (trans) {
244  (*it).pfi->real_grad_base_value(ctx, elmt_stored[k], icb != 0);
245  *mit = short_type(ctx.N());
246  }
247  else
248  elmt_stored[k] = pfp[k]->grad(ctx.ii());
249  break;
250  case GETFEM_HESSIAN_ :
251  ++mit;
252  if ((*it).pfi && (*it).pfi->target_dim() > 1) ++mit;
253  if (trans) {
254  (*it).pfi->real_hess_base_value(ctx, elmt_stored[k], icb != 0);
255  *mit = short_type(gmm::sqr(ctx.N()));
256  }
257  else {
258  base_tensor tt = pfp[k]->hess(ctx.ii());
259  aux_ind.resize(3);
260  aux_ind[2] = gmm::sqr(tt.sizes()[2]); aux_ind[1] = tt.sizes()[1];
261  aux_ind[0] = tt.sizes()[0];
262  tt.adjust_sizes(aux_ind);
263  elmt_stored[k] = tt;
264  }
265  break;
266  case GETFEM_UNIT_NORMAL_ :
267  *mit = short_type(ctx.N());
268  {
269  aux_ind.resize(1); aux_ind[0] = short_type(ctx.N());
270  elmt_stored[k].adjust_sizes(aux_ind);
271  }
272  std::copy(up.begin(), up.end(), elmt_stored[k].begin());
273  break;
274  case GETFEM_GRAD_GEOTRANS_ :
275  case GETFEM_GRAD_GEOTRANS_INV_ : {
276  size_type P = gmm::mat_ncols(ctx.K()), N=ctx.N();
277  base_matrix Bt;
278  if (it->t == GETFEM_GRAD_GEOTRANS_INV_) {
279  Bt.resize(P,N); gmm::copy(gmm::transposed(ctx.B()),Bt);
280  }
281  const base_matrix &A = (it->t==GETFEM_GRAD_GEOTRANS_) ? ctx.K():Bt;
282  aux_ind.resize(2);
283  *mit++ = aux_ind[0] = short_type(gmm::mat_nrows(A));
284  *mit = aux_ind[1] = short_type(gmm::mat_ncols(A));
285  elmt_stored[k].adjust_sizes(aux_ind);
286  std::copy(A.begin(), A.end(), elmt_stored[k].begin());
287  } break;
288  case GETFEM_NONLINEAR_ :
289  if ((*it).nl_part != 0) { /* for auxiliary fem of nonlinear_term,*/
290  /* the "prepare" method is called */
291  if ((*it).nlt->term_num() == size_type(-1)) {
292  (*it).nlt->prepare(ctx, (*it).nl_part);
293  /* the dummy assistant multiplies everybody by 1
294  -> not efficient ! */
295  }
296  aux_ind.resize(1); aux_ind[0] = 1;
297  elmt_stored[k].adjust_sizes(aux_ind); elmt_stored[k][0] = 1.;
298  } else {
299  if ((*it).nlt->term_num() == size_type(-1)) {
300  const bgeot::multi_index &nltsizes
301  = (*it).nlt->sizes(ctx.convex_num());
302  elmt_stored[k].adjust_sizes(nltsizes);
303  (*it).nlt->compute(ctx, elmt_stored[k]);
304  (*it).nlt->term_num() = k;
305  for (dim_type ii = 0; ii < nltsizes.size(); ++ii)
306  *mit++ = nltsizes[ii];
307  --mit;
308  } else {
309  elmt_stored[k] = elmt_stored[(*it).nlt->term_num()];
310  const bgeot::multi_index &nltsizes = elmt_stored[k].sizes();
311  for (dim_type ii = 0; ii < nltsizes.size(); ++ii)
312  *mit++ = nltsizes[ii];
313  --mit;
314  }
315  }
316  break;
317  }
318  }
319 
320  GMM_ASSERT1(mit == sizes.end(), "internal error");
321 
322  //expand_product_old(t,J*pai->coeff(ctx.ii()), first);
323  scalar_type c = J*pai->coeff(ctx.ii());
324  if (!icb) {
325  if (first) { t.adjust_sizes(sizes); }
326  expand_product_daxpy(t, c, first);
327  } else {
328  icb->eltm.resize(0);
329  for (unsigned k=0; k != pme->size(); ++k) {
330  if (icb && !((*pme)[k].t == GETFEM_NONLINEAR_
331  && (*pme)[k].nl_part != 0))
332  icb->eltm.push_back(&elmt_stored[k]);
333  }
334  icb->exec(t, first, c);
335  }
336  }
337 
338 
339  void expand_product_old(base_tensor &t, scalar_type J, bool first) const {
340  scalar_type V;
341  size_type k;
342  if (first) std::fill(t.begin(), t.end(), 0.0);
343  base_tensor::iterator pt = t.begin();
344  std::vector<base_tensor::const_iterator> pts(pme->size());
345  std::vector<scalar_type> Vtab(pme->size());
346  for (k = 0; k < pme->size(); ++k)
347  pts[k] = elmt_stored[k].begin();
348 
349  size_type k0 = 0;
350  unsigned n0 = unsigned(elmt_stored[0].size());
351  /*while (elmt_stored[k0].size() == 1 && k0+1 < pme->size()) {
352  J *= elmt_stored[k0][0];
353  ++k0; n0 = elmt_stored[k0].size();
354  }*/
355  base_tensor::const_iterator pts0 = pts[k0];
356 
357 
358  k = pme->size()-1; Vtab[k] = J;
359  /* very heavy expansion .. takes much time */
360  do {
361  for (V = Vtab[k]; k!=k0; --k)
362  Vtab[k-1] = V = *pts[k] * V;
363  for (k=0; k < n0; ++k)
364  *pt++ += V * pts0[k];
365  for (k=k0+1; k != pme->size() && ++pts[k] == elmt_stored[k].end(); ++k)
366  pts[k] = elmt_stored[k].begin();
367  } while (k != pme->size());
368  GMM_ASSERT1(pt == t.end(), "Internal error");
369  }
370 
371  /* do the tensorial product using the blas function daxpy (much more
372  efficient than a loop).
373 
374  efficiency is maximized when the first tensor has a large dimension
375  */
376  void expand_product_daxpy(base_tensor &t, scalar_type J, bool first)const {
377  size_type k;
378  base_tensor::iterator pt = t.begin();
379  THREAD_SAFE_STATIC std::vector<base_tensor::const_iterator> pts;
380  THREAD_SAFE_STATIC std::vector<base_tensor::const_iterator> es_beg;
381  THREAD_SAFE_STATIC std::vector<base_tensor::const_iterator> es_end;
382  THREAD_SAFE_STATIC std::vector<scalar_type> Vtab;
383 
384  pts.resize(0); pts.resize(pme->size()); // resize(0) necessary, do not remove
385  es_beg.resize(0); es_beg.resize(pme->size());
386  es_end.resize(0); es_end.resize(pme->size());
387  Vtab.resize(pme->size());
388  size_type nm = 0;
389  if (first) memset(&(*t.begin()), 0, t.size()*sizeof(*t.begin())); //std::fill(t.begin(), t.end(), 0.0);
390  for (k = 0, nm = 0; k < pme->size(); ++k) {
391  if (elmt_stored[k].size() != 1) {
392  es_beg[nm] = elmt_stored[k].begin();
393  es_end[nm] = elmt_stored[k].end();
394  pts[nm] = elmt_stored[k].begin();
395  ++nm;
396  } else J *= elmt_stored[k][0];
397  }
398  if (nm == 0) {
399  t[0] += J;
400  } else {
401  BLAS_INT n0 = BLAS_INT(es_end[0] - es_beg[0]);
402  base_tensor::const_iterator pts0 = pts[0];
403 
404  /* very heavy reduction .. takes much time */
405  k = nm-1; Vtab[k] = J;
406  BLAS_INT one = BLAS_INT(1);
407  scalar_type V;
408  do {
409  for (V = Vtab[k]; k; --k)
410  Vtab[k-1] = V = *pts[k] * V;
411  GMM_ASSERT1(pt+n0 <= t.end(), "Internal error");
412  gmm::daxpy_(&n0, &V, const_cast<double*>(&(pts0[0])), &one,
413  (double*)&(*pt), &one);
414  pt+=n0;
415  for (k=1; k != nm && ++pts[k] == es_end[k]; ++k)
416  pts[k] = es_beg[k];
417  } while (k != nm);
418  GMM_ASSERT1(pt == t.end(), "Internal error");
419  }
420  }
421 
422 
423  void pre_tensors_for_linear_trans(bool volumic) const {
424 
425  if ((volumic && volume_computed) || (!volumic && faces_computed)) return;
426  // scalar_type exectime = ftool::uclock_sec();
427 
428  bgeot::multi_index sizes = pme->sizes(0), mi(sizes.size());
429  bgeot::multi_index::iterator mit = sizes.begin(), mite = sizes.end();
430  size_type f = 1;
431  for ( ; mit != mite; ++mit, ++f) f *= *mit;
432  if (f > 1000000)
433  GMM_WARNING2("Warning, very large elementary computations.\n"
434  << "Be sure you need to compute this elementary matrix.\n"
435  << "(sizes = " << sizes << " )\n");
436 
437  base_tensor aux(sizes);
438  std::fill(aux.begin(), aux.end(), 0.0);
439  if (volumic) {
440  volume_computed = true;
441  mref[0] = aux;
442  }
443  else {
444  faces_computed = true;
445  std::fill(mref.begin()+1, mref.end(), aux);
446  }
447 
448  if (is_ppi) // pour accelerer, il faudrait pr�calculer les d�riv�es
449  {
450  base_poly P(dim, 0), Q(dim, 0), R(dim, 0);
451  size_type old_ind = size_type(-1), ind;
452  for ( ; !mi.finished(sizes); mi.incrementation(sizes)) {
453 
454  mat_elem_type::const_iterator it = pme->begin(), ite = pme->end();
455  mit = mi.begin();
456 
457  ind = *mit; ++mit;
458 
459  if ((*it).pfi) {
460  if ((*it).pfi->target_dim() > 1)
461  { ind += (*it).pfi->nb_base(0) * (*mit); ++mit; }
462 
463  Q = ((ppolyfem)((*it).pfi).get())->base()[ind];
464  }
465 
466  switch ((*it).t) {
467  case GETFEM_GRAD_ : Q.derivative(short_type(*mit)); ++mit; break;
468  case GETFEM_HESSIAN_ :
469  Q.derivative(short_type(*mit % dim));
470  Q.derivative(short_type(*mit / dim));
471  ++mit; break;
472  case GETFEM_BASE_ : break;
473  case GETFEM_GRAD_GEOTRANS_:
474  case GETFEM_GRAD_GEOTRANS_INV_:
475  case GETFEM_UNIT_NORMAL_ :
476  case GETFEM_NONLINEAR_ :
477  GMM_ASSERT1(false,
478  "Normals, gradients of geotrans and non linear "
479  "terms are not compatible with exact integration, "
480  "use an approximate method instead");
481  }
482  ++it;
483 
484  if (it != ite && *mit != old_ind) {
485  old_ind = *mit;
486  P.one();
487  for (; it != ite; ++it) {
488  ind = *mit; ++mit;
489 
490  if ((*it).pfi->target_dim() > 1)
491  { ind += (*it).pfi->nb_base(0) * (*mit); ++mit; }
492  R = ((ppolyfem)((*it).pfi).get())->base()[ind];
493 
494  switch ((*it).t) {
495  case GETFEM_GRAD_ :
496  R.derivative(short_type(*mit)); ++mit;
497  break;
498  case GETFEM_HESSIAN_ :
499  R.derivative(short_type(*mit % dim));
500  R.derivative(short_type(*mit / dim));
501  ++mit; break;
502  case GETFEM_BASE_ : break;
503  case GETFEM_UNIT_NORMAL_ :
504  case GETFEM_GRAD_GEOTRANS_:
505  case GETFEM_GRAD_GEOTRANS_INV_ :
506  case GETFEM_NONLINEAR_ :
507  GMM_ASSERT1(false, "No nonlinear term allowed here");
508  }
509  P *= R;
510  }
511  }
512  R = P * Q;
513  if (volumic) mref[0](mi) = bgeot::to_scalar(ppi->int_poly(R));
514  for (f = 0; f < nbf && !volumic; ++f)
515  mref[f+1](mi) = bgeot::to_scalar(ppi->int_poly_on_face(R, short_type(f)));
516  }
517  }
518  else {
519  bool first = true;
520  fem_interpolation_context ctx;
521  size_type ind_l = 0, nb_ptc = pai->nb_points_on_convex(),
522  nb_pt_l = nb_ptc, nb_pt_tot =(volumic ? nb_ptc : pai->nb_points());
523  for (size_type ip = (volumic ? 0:nb_ptc); ip < nb_pt_tot; ++ip) {
524  while (ip == nb_pt_l && ind_l < nbf)
525  { nb_pt_l += pai->nb_points_on_face(short_type(ind_l)); ind_l++; }
526  ctx.set_ii(ip);
527  add_elem(mref[ind_l], ctx, 1.0, first, false, NULL, sizes);
528  first = false;
529  }
530  }
531  // cout << "precompute Mat elem computation time : "
532  // << ftool::uclock_sec() - exectime << endl;
533  }
534 
535 
536  void compute(base_tensor &t, const base_matrix &G, short_type ir,
537  size_type elt, mat_elem_integration_callback *icb = 0) const {
538  dim_type P = dim_type(dim), N = dim_type(G.nrows());
539  short_type NP = short_type(pgt->nb_points());
540  fem_interpolation_context ctx(pgp, 0, 0, G, elt,
541  short_type(ir-1));
542 
543  GMM_ASSERT1(G.ncols() == NP, "dimensions mismatch");
544  if (ir > 0) {
545  up.resize(N); un.resize(P);
546  //un = pgt->normals()[ir-1];
547  gmm::copy(pgt->normals()[ir-1],un);
548  }
549  base_tensor taux;
550  bool flag = false;
551 
552  if (!computed_on_real_element) {
553  pre_tensors_for_linear_trans(ir == 0);
554  const base_matrix& B = ctx.B(); // compute B and J
555  scalar_type J=ctx.J();
556  if (ir > 0) {
557  gmm::mult(B, un, up);
558  scalar_type nup = gmm::vect_norm2(up);
559  J *= nup; //up /= nup;
560  gmm::scale(up,1.0/nup);
561  }
562 
563  t = mref[ir]; gmm::scale(t.as_vector(), J);
564 
565  if (grad_reduction.size() > 0) {
566  std::deque<short_type>::const_iterator it = grad_reduction.begin(),
567  ite = grad_reduction.end();
568  for ( ; it != ite; ++it) {
569  (flag ? t:taux).mat_transp_reduction(flag ? taux:t, B, *it);
570  flag = !flag;
571  }
572  }
573 
574  if (K_reduction.size() > 0) {
575  std::deque<short_type>::const_iterator it = K_reduction.begin(),
576  ite = K_reduction.end();
577  for ( ; it != ite; ++it) {
578  (flag ? t:taux).mat_transp_reduction(flag ? taux:t, ctx.K(), *it);
579  // (flag ? t:taux).mat_transp_reduction(flag ? taux:t, B, *it);
580  flag = !flag;
581  }
582  }
583 
584  if (hess_reduction.size() > 0) {
585  std::deque<short_type>::const_iterator it = hess_reduction.begin(),
586  ite = hess_reduction.end();
587  for (short_type l = 1; it != ite; ++it, l = short_type(l*2)) {
588  (flag ? t:taux).mat_transp_reduction(flag ? taux:t, ctx.B3(), *it);
589  flag = !flag;
590  }
591  }
592 
593  } else { // non linear transformation and methods defined on real elements
594  bgeot::multi_index sizes = pme->sizes(elt);
595 
596  bool first = true;
597  for (size_type ip=(ir == 0) ? 0 : pai->repart()[ir-1];
598  ip < pai->repart()[ir]; ++ip, first = false) {
599  ctx.set_ii(ip);
600  const base_matrix& B = ctx.B(); // J computed as side-effect
601  scalar_type J = ctx.J();
602  if (ir > 0) {
603  gmm::mult(B, un, up);
604  scalar_type nup = gmm::vect_norm2(up);
605  J *= nup; /*up /= nup;*/gmm::scale(up,1.0/nup);
606  }
607  add_elem(t, ctx, J, first, true, icb, sizes);
608  }
609 
610  // GMM_ASSERT1(!first, "No integration point on this element.");
611  if (first) {
612  GMM_WARNING3("No integration point on this element. "
613  "Caution, returning a null tensor");
614  t.adjust_sizes(sizes); gmm::clear(t.as_vector());
615  }
616  }
617 
618  /* Applying linear transformation for non tau-equivalent elements. */
619 
620  if (trans_reduction.size() > 0 && !icb) {
621  std::deque<short_type>::const_iterator it = trans_reduction.begin(),
622  ite = trans_reduction.end();
623  std::deque<pfem>::const_iterator iti = trans_reduction_pfi.begin();
624  for ( ; it != ite; ++it, ++iti) {
625  ctx.set_pf(*iti); // cout << "M = " << ctx.M() << endl;
626  (flag ? t:taux).mat_transp_reduction(flag ? taux:t, ctx.M(), *it);
627  flag = !flag;
628  }
629  }
630  if (flag) t = taux;
631  }
632 
633  void compute(base_tensor &t, const base_matrix &G, size_type elt,
634  mat_elem_integration_callback *icb) const
635  { compute(t, G, 0, elt, icb); }
636 
637  void compute_on_face(base_tensor &t, const base_matrix &G,
638  short_type f, size_type elt,
639  mat_elem_integration_callback *icb) const
640  { compute(t, G, short_type(f+1), elt, icb); }
641  };
642 
643  pmat_elem_computation mat_elem(pmat_elem_type pm, pintegration_method pi,
644  bgeot::pgeometric_trans pg,
645  bool prefer_comp_on_real_element) {
646  dal::pstatic_stored_object_key
647  pk = std::make_shared<emelem_comp_key_>(pm, pi, pg,
648  prefer_comp_on_real_element);
649  dal::pstatic_stored_object o = dal::search_stored_object(pk);
650  if (o) return std::dynamic_pointer_cast<const mat_elem_computation>(o);
651  pmat_elem_computation
652  p = std::make_shared<emelem_comp_structure_>
653  (pm, pi, pg, prefer_comp_on_real_element);
654  dal::add_stored_object(pk, p, pm, pi, pg);
655  return p;
656  }
657 
658 
659 } /* end of namespace getfem. */
660 
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.
getfem::fem_precomp
pfem_precomp fem_precomp(pfem pf, bgeot::pstored_point_tab pspt, dal::pstatic_stored_object dep)
Handles precomputations for FEM.
Definition: getfem_fem.cc:4333
getfem_mat_elem.h
elementary computations (used by the generic assembly).
getfem::ppolyfem
const typedef fem< bgeot::base_poly > * ppolyfem
Classical polynomial FEM.
Definition: getfem_fem.h:599
getfem::mat_elem
pmat_elem_computation mat_elem(pmat_elem_type pm, pintegration_method pi, bgeot::pgeometric_trans pg, bool prefer_comp_on_real_element=false)
allocate a structure for computation (integration over elements or faces of elements) of elementary t...
Definition: getfem_mat_elem.cc:643
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
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
getfem_omp.h
Tools for multithreaded, OpenMP and Boost based parallelization.
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
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::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem_fem.h
Definition of the finite element methods.