GetFEM  5.5
getfem_projected_fem.cc
1 /*===========================================================================
2 
3  Copyright (C) 2012-2026 Yves Renard, Konstantinos Poulios
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, see https://www.gnu.org/licenses/.
18 
19 ===========================================================================*/
20 
21 #include "getfem/getfem_projected_fem.h"
22 
23 namespace getfem {
24 
25  typedef bgeot::convex<base_node>::dref_convex_pt_ct dref_convex_pt_ct;
26 // typedef bgeot::basic_mesh::ref_mesh_face_pt_ct ref_mesh_face_pt_ct;
27 
28  /* calculates the projection of a point on the face of a convex
29  * Input:
30  * pgt : the geometric transformation of the convex
31  * G_cv: the nodes of the convex, stored in columns
32  * fc : the face of the convex to project on
33  * pt : the point to be projected
34  * Output:
35  * proj_ref: the projected point in the reference element
36  * proj : the projected point in the real element
37  */
38  void projection_on_convex_face
39  (const bgeot::pgeometric_trans pgt, const base_matrix &G_cv,
40  const short_type fc, const base_node &pt,
41  base_node &proj_ref, base_node &proj) {
42 
43  size_type N = gmm::mat_nrows(G_cv); // dimension of the target space
44  size_type P = pgt->dim(); // dimension of the reference element space
45 
46  size_type nb_pts_cv = gmm::mat_ncols(G_cv);
47  size_type nb_pts_fc = pgt->structure()->nb_points_of_face(fc);
48 
49  GMM_ASSERT1( N == pt.size(), "Dimensions mismatch");
50  GMM_ASSERT1( nb_pts_cv == pgt->nb_points(), "Dimensions mismatch");
51 
52  bgeot::convex_ind_ct ind_pts_fc = pgt->structure()->ind_points_of_face(fc);
53 
54  base_matrix G_fc(N, nb_pts_fc);
55  for (size_type i=0; i < nb_pts_fc; i++)
56  gmm::copy(gmm::mat_col(G_cv,ind_pts_fc[i]),gmm::mat_col(G_fc,i));
57 
58  // Local base on reference face
59  base_matrix base_ref_fc(P,P-1);
60  {
61  dref_convex_pt_ct dref_pts_fc = pgt->convex_ref()->dir_points_of_face(fc);
62  GMM_ASSERT1( dref_pts_fc.size() == P, "Dimensions mismatch");
63  for (size_type i = 0; i < P-1; ++i)
64  gmm::copy(dref_pts_fc[i+1] - dref_pts_fc[0], gmm::mat_col(base_ref_fc,i));
65  }
66 
67  proj_ref.resize(P);
68  proj.resize(N);
69 
70  base_node vres(P); // residual vector
71  scalar_type res= 1.;
72 
73  // initial guess
74  proj_ref = gmm::mean_value(pgt->convex_ref()->points_of_face(fc));
75 
76  base_vector val(nb_pts_fc);
77  pgt->poly_vector_val(proj_ref, ind_pts_fc, val);
78  gmm::mult(G_fc, val, proj);
79 
80  base_matrix K(N,P-1);
81 
82  base_matrix grad_fc(nb_pts_fc, P);
83  base_matrix grad_fc1(nb_pts_fc, P-1);
84  base_matrix B(N,P-1), BB(N,P), CS(P-1,P-1);
85 
86  scalar_type EPS = 10E-12;
87  unsigned cnt = 50;
88  while (res > EPS && --cnt) {
89  // computation of the pseudo inverse matrix B at point proj_ref
90  pgt->poly_vector_grad(proj_ref, ind_pts_fc, grad_fc);
91  gmm::mult(grad_fc, base_ref_fc, grad_fc1);
92  gmm::mult(G_fc, grad_fc1, K);
93  gmm::mult(gmm::transposed(K), K, CS);
94  gmm::lu_inverse(CS);
95  gmm::mult(K, CS, B);
96  gmm::mult(B, gmm::transposed(base_ref_fc), BB);
97 
98  // Projection onto the face of the convex
99  gmm::mult_add(gmm::transposed(BB), pt - proj, proj_ref);
100  pgt->poly_vector_val(proj_ref, ind_pts_fc, val);
101  gmm::mult(G_fc, val, proj);
102 
103  gmm::mult(gmm::transposed(BB), pt - proj, vres);
104  res = gmm::vect_norm2(vres);
105  }
106  GMM_ASSERT1( res <= EPS,
107  "Iterative pojection on convex face did not converge");
108  pgt->project_into_reference_convex(proj_ref);
109  pgt->poly_vector_val(proj_ref, ind_pts_fc, val);
110  gmm::mult(G_fc, val, proj);
111  }
112 
113 
114  /* calculates the normal at a specific point on the face of a convex
115  * Input:
116  * pgt : the geometric transformation of the convex
117  * G_cv : the nodes of the convex, stored in columns
118  * fc : the face of the convex to project on
119  * ref_pt: the point in the reference element
120  * Output:
121  * normal: the surface normal in the real element corresponding at
122  * the location of ref_pt in the reference element
123  */
124  void normal_on_convex_face
125  (const bgeot::pgeometric_trans pgt, const base_matrix &G_cv,
126  const short_type fc, const base_node &ref_pt, base_node &normal) {
127 
128  size_type N = gmm::mat_nrows(G_cv); // dimension of the target space
129  size_type P = pgt->dim(); // dimension of the reference element space
130 
131  size_type nb_pts_cv = gmm::mat_ncols(G_cv);
132  size_type nb_pts_fc = pgt->structure()->nb_points_of_face(fc);
133 
134  GMM_ASSERT1( nb_pts_cv == pgt->nb_points(), "Dimensions mismatch");
135 
136  bgeot::convex_ind_ct ind_pts_fc = pgt->structure()->ind_points_of_face(fc);
137 
138  base_matrix G_fc(N, nb_pts_fc);
139  for (size_type i=0; i < nb_pts_fc; i++)
140  gmm::copy(gmm::mat_col(G_cv,ind_pts_fc[i]),gmm::mat_col(G_fc,i));
141 
142  // Local base on reference face
143  base_matrix base_ref_fc(P,P-1);
144  {
145  dref_convex_pt_ct dref_pts_fc = pgt->convex_ref()->dir_points_of_face(fc);
146  GMM_ASSERT1( dref_pts_fc.size() == P, "Dimensions mismatch");
147  for (size_type i = 0; i < P-1; ++i)
148  gmm::copy(dref_pts_fc[i+1] - dref_pts_fc[0], gmm::mat_col(base_ref_fc,i));
149  }
150 
151  normal.resize(N);
152 
153  base_matrix K(N,P-1);
154  { // calculate K at the final point
155  base_matrix grad_fc(nb_pts_fc, P);
156  base_matrix grad_fc1(nb_pts_fc, P-1);
157  pgt->poly_vector_grad(ref_pt, ind_pts_fc, grad_fc);
158  gmm::mult(grad_fc, base_ref_fc, grad_fc1);
159  gmm::mult(G_fc, grad_fc1, K);
160  }
161 
162  base_matrix KK(N,P);
163  { // calculate KK
164  base_matrix grad_cv(nb_pts_cv, P);
165  pgt->poly_vector_grad(ref_pt, grad_cv);
166  gmm::mult(G_cv, grad_cv, KK);
167  }
168 
169  base_matrix bases_product(P-1, P);
170  gmm::mult(gmm::transposed(K), KK, bases_product);
171 
172  for (size_type i = 0; i < P; ++i) {
173  std::vector<size_type> ind(0);
174  for (size_type j = 0; j < P; ++j)
175  if (j != i ) ind.push_back(j);
176  scalar_type det = gmm::lu_det(gmm::sub_matrix(bases_product,
177  gmm::sub_interval(0, P-1),
178  gmm::sub_index(ind) ) );
179  gmm::add(gmm::scaled(gmm::mat_col(KK, i), (i % 2) ? -det : +det ), normal);
180  }
181 
182  // normalizing
183  gmm::scale(normal, 1/gmm::vect_norm2(normal));
184 
185  // ensure that normal points outwards
186  base_node cv_center(N), fc_center(N);
187  for (size_type i=0; i < nb_pts_cv; i++)
188  gmm::add(gmm::mat_col(G_cv,i), cv_center);
189  for (size_type i=0; i < nb_pts_fc; i++)
190  gmm::add(gmm::mat_col(G_fc,i), fc_center);
191  gmm::scale(cv_center, scalar_type(1)/scalar_type(nb_pts_cv));
192  gmm::scale(fc_center, scalar_type(1)/scalar_type(nb_pts_fc));
193  if (gmm::vect_sp(normal, fc_center -cv_center) < 0)
194  gmm::scale(normal, scalar_type(-1));
195  }
196 
197  /* calculates the normal at a specific point of a convex in a higher
198  * dimension space
199  * Input:
200  * pgt : the geometric transformation of the convex
201  * G_cv : the nodes of the convex, stored in columns
202  * ref_pt: the point in the reference element
203  * Output:
204  * normal: the surface normal in the real element corresponding at
205  * the location of ref_pt in the reference element
206  * (or one of the possible normals if the space dimension
207  * is more than one higher than the convex dimension)
208  */
209  void normal_on_convex
210  (const bgeot::pgeometric_trans pgt, const base_matrix &G_cv,
211  const base_node &ref_pt, base_node &normal) {
212 
213  size_type N = gmm::mat_nrows(G_cv); // dimension of the target space
214  size_type P = pgt->dim(); // dimension of the reference element space
215 
216  GMM_ASSERT1( N == 2 || N == 3, "Normal on convexes calculation is supported "
217  "only for space dimension equal to 2 or 3.");
218  GMM_ASSERT1( P < N, "Normal on convex is defined only in a space of"
219  "higher dimension.");
220 
221  size_type nb_pts = gmm::mat_ncols(G_cv);
222  base_matrix K(N,P);
223  { // calculate K at the final point
224  base_matrix grad_cv(nb_pts, P);
225  pgt->poly_vector_grad(ref_pt, grad_cv);
226  gmm::mult(G_cv, grad_cv, K);
227  }
228 
229  gmm::resize(normal,N);
230  if (P==1 && N == 2) {
231  normal[0] = -K(1,0);
232  normal[1] = K(0,0);
233  }
234  else if (P==1 && N == 3) {
235  normal[0] = K(2,0)-K(1,0);
236  normal[1] = K(0,0)-K(2,0);
237  normal[2] = K(1,0)-K(0,0);
238  }
239  else if (P==2) {
240  normal[0] = K(1,0)*K(2,1)-K(2,0)*K(1,1);
241  normal[1] = K(2,0)*K(0,1)-K(0,0)*K(2,1);
242  normal[2] = K(0,0)*K(1,1)-K(1,0)*K(0,1);
243  }
244  gmm::scale(normal, 1/gmm::vect_norm2(normal));
245  }
246 
247  void projected_fem::build_kdtree(void) const {
248  tree.clear();
249  dal::bit_vector dofs=mf_source.basic_dof_on_region(rg_source);
250  dofs.setminus(blocked_dofs);
251  dim_type qdim=target_dim();
252  for (dal::bv_visitor dof(dofs); !dof.finished(); ++dof)
253  if (dof % qdim == 0)
254  tree.add_point_with_id(mf_source.point_of_basic_dof(dof), dof);
255  }
256 
257  bool projected_fem::find_a_projected_point(const base_node &pt, base_node &ptr_proj,
258  size_type &cv_proj, short_type &fc_proj) const {
259 
260  bgeot::index_node_pair ipt;
261  //scalar_type dist =
262  tree.nearest_neighbor(ipt, pt);
263 
264  size_type cv_sel = size_type(-1);
265  short_type fc_sel = short_type(-1);
266  scalar_type dist_sel(1e10);
267  base_node proj_ref, proj_ref_sel, proj, proj_sel;
268  const getfem::mesh::ind_cv_ct cvs = mf_source.convex_to_basic_dof(ipt.i);
269  for (size_type i=0; i < cvs.size(); ++i) {
270  size_type cv = cvs[i];
271  const bgeot::pgeometric_trans pgt = mf_source.linked_mesh().trans_of_convex(cv);
272  if (rg_source.is_in(cv)) { // project on the convex
273  bool gt_invertible;
274  gic = bgeot::geotrans_inv_convex(mf_source.linked_mesh().convex(cv), pgt);
275  gic.invert(pt, proj_ref, gt_invertible);
276  if (gt_invertible) {
277  pgt->project_into_reference_convex(proj_ref);
278  proj = pgt->transform(proj_ref, mf_source.linked_mesh().points_of_convex(cv));
279  scalar_type dist = gmm::vect_dist2(pt, proj);
280  if (dist < dist_sel) {
281  dist_sel = dist;
282  cv_sel = cv;
283  fc_sel = short_type(-1);
284  proj_ref_sel = proj_ref;
285  }
286  }
287  }
288  else { // project on convex faces
289  mesh_region::face_bitset faces = rg_source.faces_of_convex(cv);
290  if (faces.count() > 0) { // this should rarely be more than one face
291  mf_source.linked_mesh().points_of_convex(cv, G);
292  short_type nbf = mf_source.linked_mesh().nb_faces_of_convex(cv);
293  for (short_type f = 0; f < nbf; ++f) {
294  if (faces.test(f)) {
295  projection_on_convex_face(pgt, G, f, pt, proj_ref, proj);
296  scalar_type dist = gmm::vect_dist2(pt, proj);
297  if (dist < dist_sel) {
298  dist_sel = dist;
299  cv_sel = cv;
300  fc_sel = f;
301  proj_ref_sel = proj_ref;
302  proj_sel = proj;
303  }
304  }
305  }
306  }
307  }
308  }
309  if (cv_sel != size_type(-1)) {
310  scalar_type elm_size = mf_source.linked_mesh().convex_radius_estimate(cv_sel);
311  if (dist_sel < 0.05*elm_size) { //FIXME
312  cv_proj = cv_sel;
313  fc_proj = fc_sel;
314  ptr_proj = proj_ref_sel;
315  return true;
316  }
317  }
318  return false;
319  }
320 
321  void projected_fem::update_from_context() const {
322  fictx_cv = size_type(-1);
323  dim_ = dim_type(-1);
324 
325  dim_type N = mf_source.linked_mesh().dim();
326  GMM_ASSERT1( N == mim_target.linked_mesh().dim(),
327  "Dimensions mismatch between the source and the target meshes");
328 
329  build_kdtree();
330 
331  elements.clear();
332  ind_dof.resize(mf_source.nb_basic_dof());
333  std::fill(ind_dof.begin(), ind_dof.end(), size_type(-1));
334  size_type max_dof = 0;
335  if (rg_target.id() != mesh_region::all_convexes().id() &&
336  rg_target.is_empty()) {
337  dim_ = mim_target.linked_mesh().dim();
338  return;
339  }
340 
341  for (mr_visitor i(rg_target); !i.finished(); ++i) {
342  size_type cv = i.cv(); // refers to the target mesh
343  short_type f = i.f(); // refers to the target mesh
344 
345  dim_type dim__ = mim_target.linked_mesh().structure_of_convex(cv)->dim();
346  if (dim_ == dim_type(-1)) {
347  dim_ = dim__;
348  if (i.is_face()) dim__ = dim_type(dim__ - 1);
349  GMM_ASSERT1(dim__ < N, "The projection should take place in lower "
350  "dimensions than the mesh dimension. Otherwise "
351  "use the interpolated_fem object instead.");
352  }
353  else
354  GMM_ASSERT1(dim_ == dim__,
355  "Convexes/faces of different dimension in the target mesh");
356 
357  pintegration_method pim = mim_target.int_method_of_element(cv);
358  GMM_ASSERT1(pim->type() == IM_APPROX,
359  "You have to use approximated integration to project a fem");
360  papprox_integration pai = pim->approx_method();
361  bgeot::pgeometric_trans pgt = mim_target.linked_mesh().trans_of_convex(cv);
362  bgeot::pgeotrans_precomp pgp =
363  bgeot::geotrans_precomp(pgt, pai->pintegration_points(), 0);
364  size_type last_cv = size_type(-1); // refers to the source mesh
365  short_type last_f = short_type(-1); // refers to the source mesh
366  size_type nb_pts = i.is_face() ? pai->nb_points_on_face(f)
367  : pai->nb_points_on_convex();
368  size_type start_pt = i.is_face() ? pai->ind_first_point_on_face(f) : 0;
369  elt_projection_data &e = elements[cv];
370  base_node gpt(N);
371  dal::bit_vector new_dofs;
372  for (size_type k = 0; k < nb_pts; ++k) {
373  pgp->transform(mim_target.linked_mesh().points_of_convex(cv),
374  start_pt + k, gpt);
375  gausspt_projection_data &gppd = e.gausspt[start_pt + k];
376  gppd.iflags = find_a_projected_point(gpt, gppd.ptref, gppd.cv, gppd.f) ? 1 : 0;
377  if (gppd.iflags) {
378  // calculate gppd.normal
379  const bgeot::pgeometric_trans pgt_source = mf_source.linked_mesh().trans_of_convex(gppd.cv);
380  mf_source.linked_mesh().points_of_convex(gppd.cv, G);
381  if (gppd.f != short_type(-1))
382  normal_on_convex_face(pgt_source, G, gppd.f, gppd.ptref, gppd.normal);
383  else
384  normal_on_convex(pgt_source, G, gppd.ptref, gppd.normal);
385  // calculate gppd.gap
386  base_node ppt = pgt_source->transform(gppd.ptref, G);
387  gppd.gap = gmm::vect_sp(gpt-ppt, gppd.normal);
388  }
389 
390  if (gppd.iflags && (last_cv != gppd.cv || last_f != gppd.f)) {
391  if (gppd.f == short_type(-1)) { // convex
392  size_type nbdof = mf_source.nb_basic_dof_of_element(gppd.cv);
393  for (size_type loc_dof = 0; loc_dof < nbdof; ++loc_dof) {
394  size_type dof = mf_source.ind_basic_dof_of_element(gppd.cv)[loc_dof];
395  if (!(blocked_dofs[dof]))
396  new_dofs.add(dof);
397  }
398  }
399  else { // convex face
400  size_type nbdof = mf_source.nb_basic_dof_of_face_of_element(gppd.cv, gppd.f);
401  for (size_type loc_dof = 0; loc_dof < nbdof; ++loc_dof) {
402  size_type dof = mf_source.ind_basic_dof_of_face_of_element(gppd.cv, gppd.f)[loc_dof];
403  if (!(blocked_dofs[dof]))
404  new_dofs.add(dof);
405  }
406  }
407  last_cv = gppd.cv;
408  last_f = gppd.f;
409  }
410  }
411 
412  size_type cnt(0);
413  dal::bit_vector old_dofs;
414  for (const size_type dof : e.inddof) {
415  old_dofs.add(dof);
416  ind_dof[dof] = cnt++;
417  }
418  for (dal::bv_visitor dof(new_dofs); !dof.finished(); ++dof)
419  if (!(old_dofs[dof])) {
420  ind_dof[dof] = cnt++;
421  e.inddof.push_back(dof);
422  }
423 
424  e.pim = pim;
425  e.nb_dof = e.inddof.size();
426  max_dof = std::max(max_dof, e.nb_dof);
427  for (size_type k = 0; k < nb_pts; ++k) {
428  gausspt_projection_data &gppd = e.gausspt[start_pt + k];
429  if (gppd.iflags) {
430  if (gppd.f == short_type(-1)) { // convex
431  size_type nbdof = mf_source.nb_basic_dof_of_element(gppd.cv);
432  for (size_type loc_dof = 0; loc_dof < nbdof; ++loc_dof) {
433  size_type dof = mf_source.ind_basic_dof_of_element(gppd.cv)[loc_dof];
434  gppd.local_dof[loc_dof] = new_dofs.is_in(dof) ? ind_dof[dof]
435  : size_type(-1);
436  }
437  }
438  else { // convex face
439  size_type nbdof = mf_source.nb_basic_dof_of_face_of_element(gppd.cv, gppd.f);
440  pfem pf = mf_source.fem_of_element(gppd.cv);
441  bgeot::convex_ind_ct ind_pts_fc = pf->structure(gppd.cv)->ind_points_of_face(gppd.f);
442  unsigned rdim = target_dim() / pf->target_dim();
443  if (rdim == 1)
444  for (size_type loc_dof = 0; loc_dof < nbdof; ++loc_dof) { // local dof with respect to the source convex face
445  size_type dof = mf_source.ind_basic_dof_of_face_of_element(gppd.cv, gppd.f)[loc_dof];
446  size_type loc_dof2 = ind_pts_fc[loc_dof]; // local dof with respect to the source convex
447  gppd.local_dof[loc_dof2] = new_dofs.is_in(dof) ? ind_dof[dof]
448  : size_type(-1);
449  }
450  else
451  for (size_type ii = 0; ii < nbdof/rdim; ++ii)
452  for (size_type jj = 0; jj < rdim; ++jj) {
453  size_type loc_dof = ii*rdim + jj; // local dof with respect to the source convex face
454  size_type dof = mf_source.ind_basic_dof_of_face_of_element(gppd.cv, gppd.f)[loc_dof];
455  size_type loc_dof2 = ind_pts_fc[ii]*rdim + jj; // local dof with respect to the source convex
456  gppd.local_dof[loc_dof2] = new_dofs.is_in(dof) ? ind_dof[dof]
457  : size_type(-1);
458  }
459  }
460  }
461  }
462 
463  for (const size_type dof : e.inddof)
464  ind_dof[dof] = size_type(-1);
465 
466  }
467  /** setup global dofs, with dummy coordinates */
468  base_node P(dim()); gmm::fill(P,1./20);
469  std::vector<base_node> node_tab_(max_dof, P);
470  pspt_override = bgeot::store_point_tab(node_tab_);
471  pspt_valid = false;
472  dof_types_.resize(max_dof);
473  std::fill(dof_types_.begin(), dof_types_.end(),
474  global_dof(dim()));
475 
476  /* ind_dof should be kept filled with -1 ( real_base_value and
477  grad_base_value expect that)
478  */
479  std::fill(ind_dof.begin(), ind_dof.end(), size_type(-1));
480  }
481 
482  size_type projected_fem::nb_dof(size_type cv) const
483  {
484  context_check();
485  GMM_ASSERT1(mim_target.linked_mesh().convex_index().is_in(cv),
486  "Wrong convex number: " << cv);
487  std::map<size_type,elt_projection_data>::const_iterator eit;
488  eit = elements.find(cv);
489  return (eit != elements.end()) ? eit->second.nb_dof : 0;
490  }
491 
492  size_type projected_fem::index_of_global_dof(size_type cv, size_type i) const
493  {
494  std::map<size_type,elt_projection_data>::const_iterator eit;
495  eit = elements.find(cv);
496  GMM_ASSERT1(eit != elements.end(), "Wrong convex number: " << cv);
497  return eit->second.inddof[i];
498  }
499 
500  bgeot::pconvex_ref projected_fem::ref_convex(size_type cv) const
501  { return mim_target.int_method_of_element(cv)->approx_method()->ref_convex(); }
502 
503  const bgeot::convex<base_node> &projected_fem::node_convex(size_type cv) const
504  {
505  GMM_ASSERT1(mim_target.linked_mesh().convex_index().is_in(cv),
506  "Wrong convex number: " << cv);
507  return *(bgeot::generic_dummy_convex_ref
508  (dim(), nb_dof(cv),
509  mim_target.linked_mesh().structure_of_convex(cv)->nb_faces()));
510  }
511 
512  void projected_fem::base_value(const base_node &, base_tensor &) const
513  { GMM_ASSERT1(false, "No base values, real only element."); }
514  void projected_fem::grad_base_value(const base_node &,
515  base_tensor &) const
516  { GMM_ASSERT1(false, "No grad values, real only element."); }
517  void projected_fem::hess_base_value(const base_node &,
518  base_tensor &) const
519  { GMM_ASSERT1(false, "No hess values, real only element."); }
520 
521  inline void projected_fem::actualize_fictx(pfem pf, size_type cv,
522  const base_node &ptr) const {
523  if (fictx_cv != cv) {
524  mf_source.linked_mesh().points_of_convex(cv, G);
525  fictx = fem_interpolation_context
526  (mf_source.linked_mesh().trans_of_convex(cv), pf, base_node(), G, cv);
527  fictx_cv = cv;
528  }
529  fictx.set_xref(ptr);
530  }
531 
532  void projected_fem::real_base_value(const fem_interpolation_context& c,
533  base_tensor &t, bool) const {
534  std::map<size_type,elt_projection_data>::iterator eit;
535  eit = elements.find(c.convex_num());
536  if (eit == elements.end()) {
537  mi2[1] = target_dim();
538  mi2[0] = short_type(0);
539  t.adjust_sizes(mi2);
540  std::fill(t.begin(), t.end(), scalar_type(0));
541  return;
542  }
543 // GMM_ASSERT1(eit != elements.end(), "Wrong convex number: " << c.convex_num());
544  elt_projection_data &e = eit->second;
545 
546  mi2[1] = target_dim();
547  mi2[0] = short_type(e.nb_dof);
548  t.adjust_sizes(mi2);
549  std::fill(t.begin(), t.end(), scalar_type(0));
550  if (e.nb_dof == 0) return;
551 
552  std::map<size_type,gausspt_projection_data>::iterator git;
553  git = e.gausspt.find(c.ii());
554  if (c.have_pgp() &&
555  (c.pgp()->get_ppoint_tab()
556  == e.pim->approx_method()->pintegration_points()) &&
557  git != e.gausspt.end()) {
558  gausspt_projection_data &gppd = git->second;
559  if (gppd.iflags & 1) {
560  if (gppd.iflags & 2) {
561  t = gppd.base_val;
562  return;
563  }
564  size_type cv = gppd.cv;
565  pfem pf = mf_source.fem_of_element(cv);
566  actualize_fictx(pf, cv, gppd.ptref);
567  pf->real_base_value(fictx, taux);
568  unsigned rdim = target_dim() / pf->target_dim();
569  std::map<size_type,size_type>::const_iterator ii;
570  if (rdim == 1) // mdim == 0
571  for (size_type i = 0; i < pf->nb_dof(cv); ++i) {
572  ii = gppd.local_dof.find(i);
573  if (ii != gppd.local_dof.end() && ii->second != size_type(-1))
574  for (size_type j = 0; j < target_dim(); ++j)
575  t(ii->second,j) = taux(i,j);
576  }
577  else // mdim == 1
578  for (size_type i = 0; i < pf->nb_dof(cv); ++i)
579  for (size_type j = 0; j < target_dim(); ++j) {
580  ii = gppd.local_dof.find(i*rdim+j);
581  if (ii != gppd.local_dof.end() && ii->second != size_type(-1))
582  t(ii->second,j) = taux(i,0);
583  }
584 
585  if (store_values) {
586  gppd.base_val = t;
587  gppd.iflags |= 2;
588  }
589  }
590  }
591  else {
592  size_type cv;
593  short_type f;
594  if (find_a_projected_point(c.xreal(), ptref, cv, f)) {
595  pfem pf = mf_source.fem_of_element(cv);
596  actualize_fictx(pf, cv, ptref);
597  pf->real_base_value(fictx, taux);
598 
599  for (size_type i = 0; i < e.nb_dof; ++i)
600  ind_dof.at(e.inddof[i]) = i;
601 
602  unsigned rdim = target_dim() / pf->target_dim();
603  if (rdim == 1) // mdim == 0
604  for (size_type i = 0; i < pf->nb_dof(cv); ++i) {
605  size_type ii = ind_dof.at(mf_source.ind_basic_dof_of_element(cv)[i]);
606  if (ii != size_type(-1)) {
607  for (size_type j = 0; j < target_dim(); ++j)
608  t(ii,j) = taux(i,j);
609  }
610  }
611  else // mdim == 1
612  for (size_type i = 0; i < pf->nb_dof(cv); ++i)
613  for (size_type j = 0; j < target_dim(); ++j) {
614  size_type ij = ind_dof.at(mf_source.ind_basic_dof_of_element(cv)[i*rdim+j]);
615  if (ij != size_type(-1))
616  t(ij,j) = taux(i,0);
617  }
618 
619  for (size_type i = 0; i < e.nb_dof; ++i)
620  ind_dof[e.inddof[i]] = size_type(-1);
621  }
622  }
623 
624  }
625 
626  void projected_fem::real_grad_base_value(const fem_interpolation_context& c,
627  base_tensor &t, bool) const {
628  std::map<size_type,elt_projection_data>::iterator eit;
629  eit = elements.find(c.convex_num());
630  if (eit == elements.end()) {
631  mi3[2] = mf_source.linked_mesh().dim();
632  mi3[1] = target_dim();
633  mi3[0] = short_type(0);
634  t.adjust_sizes(mi3);
635  std::fill(t.begin(), t.end(), scalar_type(0));
636  return;
637  }
638 // GMM_ASSERT1(eit != elements.end(), "Wrong convex number: " << c.convex_num());
639  elt_projection_data &e = eit->second;
640 
641  size_type N = mf_source.linked_mesh().dim();
642  mi3[2] = short_type(N);
643  mi3[1] = target_dim();
644  mi3[0] = short_type(e.nb_dof);
645  t.adjust_sizes(mi3);
646  std::fill(t.begin(), t.end(), scalar_type(0));
647  if (e.nb_dof == 0) return;
648 
649  std::map<size_type,gausspt_projection_data>::iterator git;
650  git = e.gausspt.find(c.ii());
651  if (c.have_pgp() &&
652  (c.pgp()->get_ppoint_tab()
653  == e.pim->approx_method()->pintegration_points()) &&
654  git != e.gausspt.end()) {
655  gausspt_projection_data &gppd = git->second;
656  if (gppd.iflags & 1) {
657  if (gppd.iflags & 4) {
658  t = gppd.grad_val;
659  return;
660  }
661  size_type cv = gppd.cv;
662  pfem pf = mf_source.fem_of_element(cv);
663  actualize_fictx(pf, cv, gppd.ptref);
664  pf->real_grad_base_value(fictx, taux);
665 
666  unsigned rdim = target_dim() / pf->target_dim();
667  std::map<size_type,size_type>::const_iterator ii;
668  if (rdim == 1) // mdim == 0
669  for (size_type i = 0; i < pf->nb_dof(cv); ++i) {
670  ii = gppd.local_dof.find(i);
671  if (ii != gppd.local_dof.end() && ii->second != size_type(-1))
672  for (size_type j = 0; j < target_dim(); ++j)
673  for (size_type k = 0; k < N; ++k)
674  t(ii->second, j, k) = taux(i, j, k);
675  }
676  else // mdim == 1
677  for (size_type i = 0; i < pf->nb_dof(cv); ++i)
678  for (size_type j = 0; j < target_dim(); ++j) {
679  ii = gppd.local_dof.find(i*rdim+j);
680  if (ii != gppd.local_dof.end() && ii->second != size_type(-1))
681  for (size_type k = 0; k < N; ++k)
682  t(ii->second, j, k) = taux(i, 0, k);
683  }
684  if (store_values) {
685  gppd.grad_val = t;
686  gppd.iflags |= 4;
687  }
688  }
689  }
690  else {
691  size_type cv;
692  short_type f;
693  if (find_a_projected_point(c.xreal(), ptref, cv, f)) {
694  pfem pf = mf_source.fem_of_element(cv);
695  actualize_fictx(pf, cv, ptref);
696  pf->real_grad_base_value(fictx, taux);
697  for (size_type i = 0; i < e.nb_dof; ++i)
698  ind_dof.at(e.inddof[i]) = i;
699 
700  unsigned rdim = target_dim() / pf->target_dim();
701  if (rdim == 1) // mdim == 0
702  for (size_type i = 0; i < pf->nb_dof(cv); ++i) {
703  size_type ii = ind_dof.at(mf_source.ind_basic_dof_of_element(cv)[i]);
704  if (ii != size_type(-1))
705  for (size_type j = 0; j < target_dim(); ++j)
706  for (size_type k = 0; k < N; ++k)
707  t(ii,j,k) = taux(i,j,k);
708  }
709  else // mdim == 1
710  for (size_type i = 0; i < pf->nb_dof(cv); ++i)
711  for (size_type j = 0; j < target_dim(); ++j) {
712  size_type ij = ind_dof.at(mf_source.ind_basic_dof_of_element(cv)[i*rdim+j]);
713  if (ij != size_type(-1))
714  for (size_type k = 0; k < N; ++k)
715  t(ij,j,k) = taux(i,0,k);
716  }
717 
718  for (size_type i = 0; i < e.nb_dof; ++i)
719  ind_dof[e.inddof[i]] = size_type(-1);
720  }
721  }
722  }
723 
724  void projected_fem::real_hess_base_value
725  (const fem_interpolation_context&, base_tensor &, bool) const
726  { GMM_ASSERT1(false, "Sorry, to be done."); }
727 
728  void projected_fem::projection_data(const fem_interpolation_context& c,
729  base_node &normal, scalar_type &gap) const {
730  std::map<size_type,elt_projection_data>::iterator eit;
731  eit = elements.find(c.convex_num());
732 
733  if (eit != elements.end()) {
734  elt_projection_data &e = eit->second;
735  if (e.nb_dof == 0) { // return undefined normal vector and huge gap
736  normal = base_node(c.N());
737  gap = 1e12;
738  return;
739  }
740  std::map<size_type,gausspt_projection_data>::iterator git;
741  git = e.gausspt.find(c.ii());
742  if (c.have_pgp() &&
743  (c.pgp()->get_ppoint_tab()
744  == e.pim->approx_method()->pintegration_points()) &&
745  git != e.gausspt.end()) {
746  gausspt_projection_data &gppd = git->second;
747  if (gppd.iflags & 1) {
748  normal = gppd.normal;
749  gap = gppd.gap;
750  }
751  else { // return undefined normal vector and huge gap
752  normal = base_node(c.N());
753  gap = 1e12;
754  }
755  return;
756  }
757  }
758 
759  // new projection
760  projection_data(c.xreal(), normal, gap);
761  }
762 
763  void projected_fem::projection_data(const base_node& pt,
764  base_node &normal, scalar_type &gap) const {
765  size_type cv;
766  short_type f;
767  if (find_a_projected_point(pt, ptref, cv, f)) {
768  const bgeot::pgeometric_trans pgt = mf_source.linked_mesh().trans_of_convex(cv);
769  mf_source.linked_mesh().points_of_convex(cv, G);
770  if (f != short_type(-1))
771  normal_on_convex_face(pgt, G, f, ptref, normal);
772  else
773  normal_on_convex(pgt, G, ptref, normal);
774  base_node ppt = pgt->transform(ptref, G);
775  gap = gmm::vect_sp(pt-ppt, normal);
776  }
777  else { // return undefined normal vector and huge gap
778  normal = base_node(pt.size());
779  gap = 1e12;
780  }
781 
782  }
783 
784  dal::bit_vector projected_fem::projected_convexes() const {
785  dal::bit_vector bv;
786  std::map<size_type,elt_projection_data>::const_iterator eit;
787  for (eit = elements.begin(); eit != elements.end(); ++eit) {
788  std::map<size_type,gausspt_projection_data>::const_iterator git;
789  for (git = eit->second.gausspt.begin(); git != eit->second.gausspt.end(); ++git) {
790  if (git->second.iflags)
791  bv.add(git->second.cv);
792  }
793  }
794  return bv;
795  }
796 
797  mesh_region projected_fem::projected_target_region() const {
798  context_check();
799  mesh_region projected_target;
800  for (mr_visitor v(rg_target); !v.finished(); ++v) {
801  pintegration_method pim = mim_target.int_method_of_element(v.cv());
802  papprox_integration pai = pim->approx_method();
803  size_type start_pt = v.is_face() ? pai->ind_first_point_on_face(v.f()) : 0;
804  size_type nb_pts = v.is_face() ? pai->nb_points_on_face(v.f())
805  : pai->nb_points_on_convex();
806  bool isProjectedOn = false;
807  for (size_type ip = 0; ip != nb_pts; ++ip) {
808  auto &proj_data = elements.at(v.cv()).gausspt[start_pt + ip];
809  if (proj_data.iflags) {
810  isProjectedOn = true;
811  break;
812  }
813  }
814  if (isProjectedOn) projected_target.add(v.cv(), v.f());
815  }
816  return projected_target;
817  }
818 
819  void projected_fem::gauss_pts_stats(unsigned &ming, unsigned &maxg,
820  scalar_type &meang) const {
821  std::vector<unsigned> v(mf_source.linked_mesh().nb_allocated_convex());
822  std::map<size_type,elt_projection_data>::const_iterator eit;
823  for (eit = elements.begin(); eit != elements.end(); ++eit) {
824  std::map<size_type,gausspt_projection_data>::const_iterator git;
825  for (git = eit->second.gausspt.begin(); git != eit->second.gausspt.end(); ++git) {
826  if (git->second.iflags)
827  v[git->second.cv]++;
828  }
829  }
830 
831  ming = 100000; maxg = 0; meang = 0;
832  unsigned cntg = 0;
833  for (dal::bv_visitor cv(mf_source.linked_mesh().convex_index());
834  !cv.finished(); ++cv) {
835  ming = std::min(ming, v[cv]);
836  maxg = std::max(maxg, v[cv]);
837  meang += v[cv];
838  if (v[cv] > 0) ++cntg;
839  }
840  meang /= scalar_type(cntg);
841  }
842 
843  size_type projected_fem::memsize() const {
844  size_type sz = 0;
845  sz += blocked_dofs.memsize();
846  sz += sizeof(*this);
847  sz += elements.size() * sizeof(elt_projection_data); // Wrong for std::map
848  std::map<size_type,elt_projection_data>::const_iterator eit;
849  for (eit = elements.begin(); eit != elements.end(); ++eit) {
850  sz += eit->second.gausspt.size() * sizeof(gausspt_projection_data); // Wrong for std::map
851  sz += eit->second.inddof.capacity() * sizeof(size_type);
852  std::map<size_type,gausspt_projection_data>::const_iterator git;
853  for (git = eit->second.gausspt.begin(); git != eit->second.gausspt.end(); ++git) {
854  sz += git->second.local_dof.size() * sizeof(size_type); // Wrong for std::map
855  }
856  }
857  return sz;
858  }
859 
860  projected_fem::projected_fem(const mesh_fem &mf_source_,
861  const mesh_im &mim_target_,
862  size_type rg_source_,
863  size_type rg_target_,
864  dal::bit_vector blocked_dofs_, bool store_val)
865  : mf_source(mf_source_), mim_target(mim_target_),
866  rg_source(mf_source.linked_mesh().region(rg_source_)),
867  rg_target(mim_target.linked_mesh().region(rg_target_)),
868  store_values(store_val), blocked_dofs(blocked_dofs_), mi2(2), mi3(3) {
869  this->add_dependency(mf_source);
870  this->add_dependency(mim_target);
871  is_pol = is_lag = is_standard_fem = false; es_degree = 5;
872  is_equiv = real_element_defined = true;
873  ntarget_dim = mf_source.get_qdim();
874 
875  update_from_context();
876  }
877 
878  DAL_SIMPLE_KEY(special_projfem_key, pfem);
879 
880  pfem new_projected_fem(const mesh_fem &mf_source_, const mesh_im &mim_target_,
881  size_type rg_source_, size_type rg_target_,
882  dal::bit_vector blocked_dofs_, bool store_val) {
883  pfem pf = std::make_shared<projected_fem>
884  (mf_source_, mim_target_, rg_source_, rg_target_,blocked_dofs_,store_val);
885  dal::pstatic_stored_object_key
886  pk = std::make_shared<special_projfem_key>(pf);
887  dal::add_stored_object(pk, pf);
888  return pf;
889  }
890 
891 
892 } /* end of namespace getfem. */
893 
bgeot::convex::dir_points_of_face
dref_convex_pt_ct dir_points_of_face(short_type i) const
Direct points for a given face.
Definition: bgeot_convex.h:84
bgeot::geotrans_interpolation_context::xreal
const base_node & xreal() const
coordinates of the current point, in the real convex.
Definition: bgeot_geometric_trans.cc:298
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:457
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:63
bgeot::geotrans_inv_convex::invert
bool invert(const base_node &n, base_node &n_ref, scalar_type IN_EPS=1e-12, bool project_into_element=false)
given the node on the real element, returns the node on the reference element (even if it is outside ...
Definition: bgeot_geotrans_inv.cc:60
bgeot::kdtree::add_point_with_id
void add_point_with_id(const base_node &n, size_type i)
insert a new point, with an associated number.
Definition: bgeot_kdtree.h:119
bgeot::kdtree::clear
void clear()
reset the tree, remove all points
Definition: bgeot_kdtree.h:112
bgeot::mesh_structure::convex_index
const dal::bit_vector & convex_index() const
Return the list of valid convex IDs.
Definition: bgeot_mesh_structure.h:90
bgeot::mesh_structure::nb_faces_of_convex
short_type nb_faces_of_convex(size_type ic) const
Return the number of faces of convex ic.
Definition: bgeot_mesh_structure.h:117
bgeot::mesh_structure::structure_of_convex
pconvex_structure structure_of_convex(size_type ic) const
Return the pconvex_structure of the convex ic.
Definition: bgeot_mesh_structure.h:111
bgeot::mesh_structure::nb_allocated_convex
size_type nb_allocated_convex() const
The number of convex indexes from 0 to the index of the last convex.
Definition: bgeot_mesh_structure.h:97
getfem::context_dependencies::context_check
bool context_check() const
return true if update_from_context was called
Definition: getfem_context.h:125
getfem::fem_interpolation_context
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:749
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:151
getfem::mesh_fem::ind_basic_dof_of_element
virtual ind_dof_ct ind_basic_dof_of_element(size_type cv) const
Give an array of the dof numbers a of convex.
Definition: getfem_mesh_fem.h:453
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:315
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:282
getfem::mesh_fem::nb_basic_dof
virtual size_type nb_basic_dof() const
Return the total number of basic degrees of freedom (before the optional reduction).
Definition: getfem_mesh_fem.h:559
getfem::mesh_fem::basic_dof_on_region
virtual dal::bit_vector basic_dof_on_region(const mesh_region &b) const
Get a list of basic dof lying on a given mesh_region.
Definition: getfem_mesh_fem.cc:76
getfem::mesh_fem::ind_basic_dof_of_face_of_element
virtual ind_dof_face_ct ind_basic_dof_of_face_of_element(size_type cv, short_type f) const
Give an array of the dof numbers lying of a convex face (all degrees of freedom whose associated base...
Definition: getfem_mesh_fem.h:471
getfem::mesh_fem::point_of_basic_dof
virtual base_node point_of_basic_dof(size_type cv, size_type i) const
Return the geometrical location of a degree of freedom.
Definition: getfem_mesh_fem.cc:222
getfem::mesh_fem::nb_basic_dof_of_face_of_element
virtual size_type nb_basic_dof_of_face_of_element(size_type cv, short_type f) const
Return the number of dof lying on the given convex face.
Definition: getfem_mesh_fem.h:484
getfem::mesh_fem::fem_of_element
virtual pfem fem_of_element(size_type cv) const
Return the basic fem associated with an element (if no fem is associated, the function will crash!...
Definition: getfem_mesh_fem.h:447
getfem::mesh_fem::nb_basic_dof_of_element
virtual size_type nb_basic_dof_of_element(size_type cv) const
Return the number of degrees of freedom attached to a given convex.
Definition: getfem_mesh_fem.h:497
getfem::mesh_fem::convex_to_basic_dof
virtual const mesh::ind_cv_ct & convex_to_basic_dof(size_type d) const
Return the list of convexes attached to the specified dof.
Definition: getfem_mesh_fem.cc:266
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:46
getfem::mesh_im::int_method_of_element
virtual pintegration_method int_method_of_element(size_type cv) const
return the integration method associated with an element (in no integration is associated,...
Definition: getfem_mesh_im.h:116
getfem::mesh_im::linked_mesh
const mesh & linked_mesh() const
Give a reference to the linked mesh of type mesh.
Definition: getfem_mesh_im.h:78
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:236
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:54
getfem::mesh_region::all_convexes
static mesh_region all_convexes()
provide a default value for the mesh_region parameters of assembly procedures etc.
Definition: getfem_mesh_region.h:141
getfem::mesh::convex_radius_estimate
virtual scalar_type convex_radius_estimate(size_type ic) const
Return an estimate of the convex largest dimension.
Definition: getfem_mesh.cc:460
getfem::mesh::convex
ref_convex convex(size_type ic) const
return a bgeot::convex object for the convex number ic.
Definition: getfem_mesh.h:188
getfem::projected_fem::projected_target_region
mesh_region projected_target_region() const
faces and convexes from the target region that contain at least one Gauss point that is projected by ...
Definition: getfem_projected_fem.cc:797
getfem::projected_fem::hess_base_value
void hess_base_value(const base_node &, base_tensor &) const
Give the value of all hessians (on ref.
Definition: getfem_projected_fem.cc:517
getfem::projected_fem::nb_dof
virtual size_type nb_dof(size_type cv) const
Number of degrees of freedom.
Definition: getfem_projected_fem.cc:482
getfem::projected_fem::real_base_value
void real_base_value(const fem_interpolation_context &c, base_tensor &t, bool=true) const
Give the value of all components of the base functions at the current point of the fem_interpolation_...
Definition: getfem_projected_fem.cc:532
getfem::projected_fem::projected_convexes
dal::bit_vector projected_convexes() const
return the list of convexes of the projected mesh_fem which contain at least one gauss point (should ...
Definition: getfem_projected_fem.cc:784
getfem::projected_fem::gauss_pts_stats
void gauss_pts_stats(unsigned &ming, unsigned &maxg, scalar_type &meang) const
return the min/max/mean number of gauss points in the convexes of the projected mesh_fem
Definition: getfem_projected_fem.cc:819
getfem::projected_fem::grad_base_value
void grad_base_value(const base_node &, base_tensor &) const
Give the value of all gradients (on ref.
Definition: getfem_projected_fem.cc:514
getfem::projected_fem::ref_convex
virtual bgeot::pconvex_ref ref_convex(size_type cv) const
Return the convex of the reference element.
Definition: getfem_projected_fem.cc:500
getfem::projected_fem::real_hess_base_value
void real_hess_base_value(const fem_interpolation_context &, base_tensor &, bool=true) const
Give the hessian of all components of the base functions at the current point of the fem_interpolatio...
Definition: getfem_projected_fem.cc:725
getfem::projected_fem::node_convex
virtual const bgeot::convex< base_node > & node_convex(size_type cv) const
Gives the convex representing the nodes on the reference element.
Definition: getfem_projected_fem.cc:503
getfem::projected_fem::update_from_context
virtual void update_from_context(void) const
this function has to be defined and should update the object when the context is modified.
Definition: getfem_projected_fem.cc:321
getfem::projected_fem::base_value
void base_value(const base_node &, base_tensor &) const
Give the value of all components of the base functions at the point x of the reference element.
Definition: getfem_projected_fem.cc:512
getfem::projected_fem::real_grad_base_value
void real_grad_base_value(const fem_interpolation_context &c, base_tensor &t, bool=true) const
Give the gradient of all components of the base functions at the current point of the fem_interpolati...
Definition: getfem_projected_fem.cc:626
gmm::tab_ref_index_ref
indexed array reference (given a container X, and a set of indexes I, this class provides a pseudo-co...
Definition: gmm_ref.h:303
getfem_projected_fem.h
FEM which projects a mesh_fem on a different mesh.
gmm::mult_add
void mult_add(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1790
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:556
gmm::vect_dist2
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_dist2(const V1 &v1, const V2 &v2)
Euclidean distance between two vectors.
Definition: gmm_blas.h:596
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:209
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:263
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1275
gmm::lu_det
linalg_traits< DenseMatrixLU >::value_type lu_det(const DenseMatrixLU &LU, const Pvector &pvector)
Compute the matrix determinant (via a LU factorization)
Definition: gmm_dense_lu.h:240
gmm::lu_inverse
void lu_inverse(const DenseMatrixLU &LU, const Pvector &pvector, const DenseMatrix &AInv_)
Given an LU factored matrix, build the inverse of the matrix.
Definition: gmm_dense_lu.h:211
getfem::global_dof
pdof_description global_dof(dim_type d)
Description of a global dof, i.e.
Definition: getfem_fem.cc:542
getfem::virtual_fem::target_dim
dim_type target_dim() const
dimension of the target space.
Definition: getfem_fem.h:313
getfem::fem_interpolation_context::convex_num
size_type convex_num() const
get the current convex number
Definition: getfem_fem.cc:52
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:243
getfem::virtual_fem::dim
dim_type dim() const
dimension of the reference element.
Definition: getfem_fem.h:310
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::generic_dummy_convex_ref
pconvex_ref generic_dummy_convex_ref(dim_type nc, size_type n, short_type nf)
generic convex with n global nodes
Definition: bgeot_convex_ref.cc:857
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:185
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:283
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::new_projected_fem
pfem new_projected_fem(const mesh_fem &mf_source, const mesh_im &mim_target, size_type rg_source_=size_type(-1), size_type rg_target_=size_type(-1), dal::bit_vector blocked_dofs=dal::bit_vector(), bool store_val=true)
create a new projected FEM.
Definition: getfem_projected_fem.cc:880
bgeot::index_node_pair
store a point and the associated index for the kdtree.
Definition: bgeot_kdtree.h:58