GetFEM  5.4.2
getfem_mesher.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-2020 Julien Pommier, 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 #include "getfem/getfem_mesher.h"
23 
24 namespace getfem {
25 
26  void mesher_level_set::init_grad(void) const {
27  gradient.resize(base.dim());
28  for (dim_type d=0; d < base.dim(); ++d) {
29  gradient[d] = base; gradient[d].derivative(d);
30  }
31  initialized = 1;
32  }
33 
34  void mesher_level_set::init_hess(void) const {
35  if (initialized < 1) init_grad();
36  hessian.resize(base.dim()*base.dim());
37  for (dim_type d=0; d < base.dim(); ++d) {
38  for (dim_type e=0; e < base.dim(); ++e) {
39  hessian[d*base.dim()+e] = gradient[d];
40  hessian[d*base.dim()+e].derivative(e);
41  }
42  }
43  initialized = 2;
44  }
45 
46  scalar_type mesher_level_set::grad(const base_node &P,
47  base_small_vector &G) const {
48  if (initialized < 1) init_grad();
49  gmm::resize(G, P.size());
50  for (size_type i = 0; i < P.size(); ++i)
51  G[i] = bgeot::to_scalar(gradient[i].eval(P.begin()));
52  return (*this)(P);
53  }
54 
55  void mesher_level_set::hess(const base_node &P, base_matrix &H) const {
56  if (initialized < 2) init_hess();
57  gmm::resize(H, P.size(), P.size());
58  for (size_type i = 0; i < base.dim(); ++i)
59  for (size_type j = 0; j < base.dim(); ++j) {
60  H(i,j) = bgeot::to_scalar(hessian[i*P.size()+j].eval(P.begin()));
61  }
62  }
63 
64 
65  //
66  // Exported functions
67  //
68 
69  bool try_projection(const mesher_signed_distance& dist, base_node &X,
70  bool on_surface) {
71  base_small_vector G; base_node Y = X;
72  scalar_type d = dist.grad(X, G), dmin = gmm::abs(d);
73  size_type iter(0), count_falt(0);
74  if (on_surface || d > 0.0)
75  while (iter == 0 || dmin > 1e-15 || gmm::vect_dist2(X, Y) > 1e-15 ) {
76  gmm::copy(X, Y);
77  if (++iter > 1000) {
78  GMM_WARNING4("Try projection failed, 1000 iterations\n\n");
79  return false;
80  // return (gmm::abs(d) < 1E-10); // is there a possibility to detect
81  } // the impossibility without making 1000 iterations ?
82  gmm::scale(G, -d / std::max(1E-8, gmm::vect_norm2_sqr(G)));
83  gmm::add(G, X);
84 // scalar_type alpha(1);
85 // scalar_type dn = dist(X);
86 // while (0 && gmm::abs(dn) > 2. * gmm::abs(d) && alpha > 0.1) {
87 // alpha /= 2;
88 // gmm::add(gmm::scaled(G, -alpha), X);
89 // dn = dist(X);
90 // }
91 // if (gmm::abs(alpha) < 1E-15) return (gmm::abs(d) < 1E-10);
92  d = dist.grad(X, G);
93  if (gmm::abs(d) >= dmin*0.95) ++count_falt;
94  else { count_falt = 0; dmin = gmm::abs(d); }
95  // if (count_falt > 20) return (gmm::abs(d) < 1E-10);
96  if (count_falt > 20) return false;
97  }
98  return true;
99  }
100 
101  // Try to find an intersection of a set of signed distance d_i.
102  // Newton method on v solution to d_i(X+Gv) = 0, where the column of
103  // G are the gradients of d_i.
104  bool pure_multi_constraint_projection
105  (const std::vector<const mesher_signed_distance*> &list_constraints,
106  base_node &X, const dal::bit_vector &cts) {
107  size_type nbco = cts.card(), i(0), info, N = X.size();
108  if (!nbco) return true;
109  // cout << "nbco = " << nbco << endl;
110  std::vector<const mesher_signed_distance*> ls(nbco);
111  std::vector<scalar_type> d(nbco), v(nbco);
112  gmm::lapack_ipvt ipvt(nbco);
113  gmm::col_matrix<base_node> G(N, nbco);
114  gmm::dense_matrix<scalar_type> H(nbco, nbco);
115  base_small_vector dd(N);
116  for (dal::bv_visitor ic(cts); !ic.finished(); ++ic, ++i)
117  { ls[i] = list_constraints[ic]; d[i] = -(ls[i]->grad(X, G[i])); }
118  base_node oldX, aux(N);
119  size_type iter = 0;
120  scalar_type residual(0), alpha(0), det(0);
121 
122  if (nbco == 1) {
123  // cout << "one constraint" << endl;
124  try_projection(*(ls[0]), X, true);
125  d[0] = -(ls[0]->grad(X, G[0]));
126  } else {
127  do {
128  oldX = X;
129  det = scalar_type(0); info = 1;
130  alpha = -scalar_type(1);
131 
132  if (nbco <= N) {
133  if (nbco < N)
134  gmm::mult(gmm::transposed(G), G, H);
135  else
136  gmm::copy(gmm::transposed(G), H);
137  // cout << "H = " << H << endl;
138  info = gmm::lu_factor(H, ipvt);
139  det = scalar_type(1);
140  for (i = 0; i < nbco; ++i) det *= H(i,i);
141  }
142  // cout << "G = " << G << endl;
143 
144  if (info || gmm::abs(det) < 1E-20) {
145  dal::bit_vector cts_red = cts;
146  int eliminated = 0;
147  i = 0;
148  for (dal::bv_visitor ic(cts); !ic.finished(); ++ic, ++i) {
149  for (size_type j = 0; j < i; ++j)
150  gmm::add(gmm::scaled(G[j], -gmm::vect_sp(G[j], G[i])), G[i]);
151  scalar_type norm_gi = gmm::vect_norm2(G[i]);
152  // cout << "norm_gi = " << norm_gi << endl;
153  if (norm_gi < 1E-10)
154  { cts_red[ic] = false; eliminated++; }
155  else
156  gmm::scale(G[i], scalar_type(1)/norm_gi);
157  }
158  // cout << "G after = " << G << endl;
159  if (eliminated >= 1) {
160  // cout << "rec call with " << eliminated << " eliminated constraints" << endl;
161  pure_multi_constraint_projection(list_constraints, X, cts_red);
162  for (i = 0; i < nbco; ++i) d[i] = -(ls[i]->grad(X, G[i]));
163 
164  det = scalar_type(0); info = 1;
165  if (nbco <= N) {
166  if (nbco < N)
167  gmm::mult(gmm::transposed(G), G, H);
168  else
169  gmm::copy(G, H);
170  info = gmm::lu_factor(H, ipvt);
171  det = scalar_type(1);
172  for (i = 0; i < nbco; ++i) det *= H(i,i);
173  }
174 
175  }
176  }
177 
178  if (gmm::vect_norm2(d) > 1e-11) {
179  if (info || gmm::abs(det) < 1E-20) {
180  for (i = 0; i < nbco; ++i)
181  try_projection(*(ls[i]), X, true);
182  for (i = 0; i < nbco; ++i) d[i] = -(ls[i]->grad(X, G[i]));
183  }
184  else {
185  gmm::lu_solve(H, ipvt, v, d);
186  if (nbco < N)
187  gmm::mult(G, v, dd);
188  else
189  gmm::copy(v, dd);
190  GMM_ASSERT1(nbco <= N, "internal error");
191  gmm::add(dd, X);
192  for (i = 0; i < nbco; ++i) d[i] = -(ls[i]->grad(X, G[i]));
193  alpha = scalar_type(1);
194  if (iter > 0)
195  while (gmm::vect_norm2(d) > residual && alpha > 1E-10) {
196  alpha /= scalar_type(2);
197  gmm::add(gmm::scaled(dd, -alpha), X);
198  for (i = 0; i < nbco; ++i) d[i] = -(ls[i]->grad(X, G[i]));
199  }
200  if (alpha < 1E-15) break;
201  }
202  }
203 
204  ++iter;
205  residual = gmm::vect_norm2(d);
206  // cout << "residual = " << residual << " alpha = " << alpha;
207  // cout << " gmm::vect_dist2(oldX,X) = " << gmm::vect_dist2(oldX,X) << endl;
208  } while (residual > 1e-14 && (residual > 1e-11 || iter < 15)
209  && iter < 200);
210  // cout << "final residual = " << residual << endl;
211  }
212  for (i = 0; i < nbco; ++i) if (gmm::abs(d[i]) > SEPS) {
213  //cout << "PURE MULTI HAS FAILED for " << cts << " nb iter = " << iter << endl;
214  return false;
215  }
216  return true;
217  }
218 
219 
220 
221 
222  // return the eigenvalue of maximal abolute value for a
223  // symmetric base_matrix
224  static scalar_type max_vp(const base_matrix& M) {
225  size_type m = mat_nrows(M);
226  GMM_ASSERT1(is_hermitian(M), "Matrix is not symmetric");
227  std::vector<scalar_type> eig(m);
228  gmm::symmetric_qr_algorithm(M, eig);
229  scalar_type emax(0);
230  for (size_type i = 0; i < m; ++i) emax = std::max(emax, gmm::abs(eig[i]));
231  return emax;
232  }
233 
234  scalar_type curvature_radius_estimate(const mesher_signed_distance &dist,
235  base_node X, bool proj) {
236  if (proj) try_projection(dist, X, true);
237  base_small_vector V;
238  base_matrix H;
239  dist.grad(X, V);
240  dist.hess(X, H);
241  return gmm::vect_norm2(V) / std::max(1E-10, max_vp(H));
242  }
243 
244  scalar_type min_curvature_radius_estimate
245  (const std::vector<const mesher_signed_distance*> &list_constraints,
246  const base_node &X, const dal::bit_vector &cts, size_type hide_first) {
247  scalar_type r0 = 1E+10;
248  for (dal::bv_visitor j(cts); !j.finished(); ++j)
249  if (j >= hide_first) {
250  scalar_type r
251  = curvature_radius_estimate(*(list_constraints[j]), X);
252  r0 = std::min(r, r0);
253  }
254  return r0;
255  }
256 
257 
258  //
259  // local functions
260  //
261 
262 
263  template <typename MAT, typename MAT2> void
264  Frobenius_condition_number_sqr_gradient(const MAT& M, MAT2& G) {
265  typedef typename gmm::linalg_traits<MAT>::value_type T;
266  typedef typename gmm::number_traits<T>::magnitude_type R;
267 
268  size_type n = mat_ncols(M);
269  gmm::dense_matrix<T> B(n,n), C(n,n);
270  gmm::mult(gmm::transposed(M), M, B);
271  R trB = gmm::mat_trace(B);
272  bgeot::lu_inverse(&(*(B.begin())), n);
273  R trBinv = gmm::mat_trace(B);
274  gmm::mult(B,B,C);
275  gmm::mult(gmm::scaled(M, T(-2)*trB), C, G);
276  gmm::add(gmm::scaled(M, T(2)*trBinv), G);
277  }
278 
279 
280  struct pt_attribute {
281  bool fixed;
282  dal::bit_vector constraints;
283  bool operator<(const pt_attribute &other) const {
284  if (fixed && !other.fixed) return true;
285  else if (!fixed && other.fixed) return false;
286  else {
287  if (constraints.last_true() > other.constraints.last_true())
288  return false;
289  else if (constraints.last_true() < other.constraints.last_true())
290  return true;
291  else if (constraints.card() > other.constraints.card()) return true;
292  else if (constraints.card() < other.constraints.card()) return false;
293  else for (dal::bv_visitor i1(constraints), i2(other.constraints);
294  !i1.finished(); ++i1, ++i2) {
295  if (i1 < i2) return true;
296  else if (i2 > i1) return false;
297  }
298  }
299  return false;
300  }
301  };
302 
303  struct mesher {
304  pmesher_signed_distance dist;
305  const mesher_virtual_function& edge_len;
306  scalar_type h0, dist_point_hull, boundary_threshold_flatness;
307  size_type N, K, iter_max, iter_wtcc;
308  int prefind, noisy;
309  base_node bounding_box_min, bounding_box_max;
310  base_vector L, L0;
311 
312  std::vector<base_node> pts, pts_prev;
313  std::vector<const pt_attribute*> pts_attr;
314  std::set<pt_attribute> attributes_set;
315  gmm::dense_matrix<size_type> t;
316 
317  scalar_type ptol, ttol, L0mult, deltat, geps, deps;
318 
319  std::vector<const mesher_signed_distance*> constraints;
320 
321  gmm::dense_matrix<scalar_type> W;
322  bgeot::mesh_structure edges_mesh;
323 
324  std::vector<size_type> attracted_points;
325  std::vector<base_node> attractor_points;
326 
327  mesher(size_type K_,
328  const pmesher_signed_distance &dist_,
329  const mesher_virtual_function &edge_len_,
330  scalar_type h0_,
331  mesh &m, const std::vector<base_node> &fixed_points,
332  int noise, size_type itm, int pref,
333  scalar_type dph, scalar_type btf)
334  : dist(dist_), edge_len(edge_len_), dist_point_hull(dph),
335  boundary_threshold_flatness(btf), iter_max(itm), prefind(pref),
336  noisy(noise) {
337  if (noise == -1) noisy = gmm::traces_level::level() - 2;
338  K=K_; h0=h0_;
339  ptol = 0.0025;
340  ttol = .1;
341  dist->bounding_box(bounding_box_min,bounding_box_max);
342  N = bounding_box_min.size();
343  if (N == 2) {
344  L0mult = 1.2; deltat = .2; geps = .001*h0;
345  } else {
346  L0mult=1+0.4/pow(scalar_type(2),scalar_type(N-1));
347  deltat=.1; geps=1e-1*h0;
348  }
349  deps=sqrt(1e-8)*h0;
350  dist->register_constraints(this->constraints);
351 
352  bgeot::pgeometric_trans pgt = bgeot::simplex_geotrans(N,1);
353  gmm::resize(W,N,N);
354  base_matrix G(N,N+1);
355  vectors_to_base_matrix(G,
356  bgeot::equilateral_simplex_of_reference(dim_type(N))->points());
357  gmm::mult(G, bgeot::geotrans_precomp(pgt, pgt->convex_ref()->pspt(),
358  0)->grad(0), W);
359  bgeot::lu_inverse(&(*(W.begin())), N);
360  do_build_mesh(m,fixed_points);
361  }
362 
363  template<class TAB> scalar_type simplex_quality(const TAB &ppts) {
364  base_matrix G(N,N), GW(N, N);
365  for (size_type i=0; i < N; ++i) {
366  base_node P = ppts[i+1] - ppts[0];
367  std::copy(P.const_begin(), P.const_begin()+N, G.begin()+i*N);
368  }
369  gmm::mult(G, W, GW);
370  return gmm::abs(1./gmm::condition_number(GW));
371  }
372 
373  scalar_type worst_element, best_element;
374 
375  scalar_type fbcond_cost_function(const base_vector &c) {
376  unsigned nbt = unsigned(gmm::mat_ncols(t));
377  scalar_type cost = 0;
378  base_matrix S(N,N), SW(N,N);
379  worst_element = 1.; best_element = 1E40;
380  for (unsigned i=0; i < nbt; ++i) {
381  for (size_type j=0; j < N; ++j) {
382  for (size_type k=0; k < N; ++k) {
383  S(k,j) = c[t(j+1,i)*N+k] - c[t(0,i)*N+k];
384  }
385  }
386  gmm::mult(S,W,SW);
387  if (bgeot::lu_det(&(*(SW.begin())), N) < 1E-16) cost += 1e30;
388  else {
389  scalar_type qual = gmm::Frobenius_condition_number_sqr(SW);
390  cost += qual;
391  worst_element = std::max(worst_element, qual / scalar_type(N*N));
392  best_element = std::min(best_element, qual / scalar_type(N*N));
393  }
394  }
395  return cost / scalar_type(N * N);
396  }
397 
398  void fbcond_cost_function_derivative(const base_vector& c,
399  base_vector &grad) {
400  gmm::clear(grad);
401  base_matrix Dcond(N,N), G(N,N), S(N,N), SW(N,N);
402  unsigned nbt = unsigned(gmm::mat_ncols(t));
403 
404  for (unsigned i=0; i < nbt; ++i) {
405  for (size_type j=0; j < N; ++j) {
406  for (size_type k=0; k < N; ++k) {
407  S(k,j) = c[t(j+1,i)*N+k] - c[t(0,i)*N+k];
408  }
409  }
410  gmm::mult(S,W,SW);
411  Frobenius_condition_number_sqr_gradient(SW,Dcond);
412  gmm::mult(Dcond, gmm::transposed(W), G);
413  // gmm::mult(W, gmm::transposed(Dcond), G);
414  for (size_type j=0; j < N; ++j) {
415  for (size_type k=0; k < N; ++k) {
416  grad[t(j+1,i)*N+k] += G(k,j);
417  grad[t(0,i)*N+k] -= G(k,j);
418  }
419  }
420  }
421  for (unsigned i=0; i < pts.size(); ++i) {
422  if (pts_attr[i]->constraints.card() || pts_attr[i]->fixed)
423  for (size_type k=0; k < N; ++k) {
424  grad[i*N+k] = 0;
425  }
426  }
427  gmm::scale(grad, scalar_type(1) / scalar_type(N * N));
428 
429  }
430 
431  struct fbcond_cost_function_object {
432  mesher &m;
433  fbcond_cost_function_object(mesher &m_) : m(m_) {}
434  scalar_type operator()(const base_vector& c) const
435  { return m.fbcond_cost_function(c); }
436  };
437 
438  struct fbcond_cost_function_derivative_object {
439  mesher &m;
440  fbcond_cost_function_derivative_object(mesher &m_) : m(m_) {}
441  void operator()(const base_vector& c, base_vector &grad) const
442  { m.fbcond_cost_function_derivative(c, grad); }
443  };
444 
445  void optimize_quality() {
446 
447  size_type nbt = gmm::mat_ncols(t);
448 
449  if (noisy > 0) cout << "Quality post-optimization\n";
450  base_vector X(pts.size() * N);
451  for (unsigned i=0; i < pts.size(); ++i)
452  gmm::copy_n(pts[i].const_begin(), N, X.begin() + i*N);
453 
454  base_matrix S(N,N), SW(N,N);
455  for (unsigned i=0; i < nbt; ++i) {
456  for (size_type j=0; j < N; ++j) {
457  for (size_type k=0; k < N; ++k) {
458  S(k,j) = X[t(j+1,i)*N+k] - X[t(0,i)*N+k];
459  }
460  }
461  if (bgeot::lu_det(&(*(S.begin())), N) < 0) {
462  std::swap(t(0,i), t(1,i));
463  for (size_type j=0; j < N; ++j) {
464  for (size_type k=0; k < N; ++k) {
465  S(k,j) = X[t(j+1,i)*N+k] - X[t(0,i)*N+k];
466  }
467  }
468  }
469  if (noisy > 0 && gmm::abs(bgeot::lu_det(&(*(S.begin())), N)) < 1e-10)
470  cout << "Element " << i << " is very bad, det = "
471  << gmm::abs(lu_det(S)) << "\n";
472  gmm::mult(S,W,SW);
473  }
474 
475  if (noisy > 0) {
476  cout << "Initial quality: "
477  << fbcond_cost_function(X)/scalar_type(nbt);
478  cout << ", best element: " << sqrt(best_element)
479  << " worst element: " << sqrt(worst_element) << endl;
480  }
481  gmm::iteration iter; iter.set_noisy(noisy-1); iter.set_maxiter(1000);
482  iter.set_resmax(5E-2*sqrt(scalar_type(nbt)));
483  gmm::bfgs(fbcond_cost_function_object(*this),
484  fbcond_cost_function_derivative_object(*this),
485  X, 10, iter, 0, 0.001, float(gmm::mat_ncols(t)));
486 
487  if (noisy > 0) cout << "Final quality: "
488  << fbcond_cost_function(X)/scalar_type(nbt)
489  << ", best element: " << sqrt(best_element)
490  << " worst element: " << sqrt(worst_element) << endl;
491 
492  for (unsigned i=0; i < pts.size(); ++i)
493  gmm::copy_n(X.begin() + i*N, N, pts[i].begin());
494  }
495 
496  void delete_element(size_type ic) {
497  if (ic != gmm::mat_ncols(t)-1) {
498  for (size_type k=0; k < N+1; ++k)
499  std::swap(t(k,ic), t(k,gmm::mat_ncols(t)-1));
500  }
501  t.resize(N+1,gmm::mat_ncols(t)-1);
502  }
503 
504  scalar_type quality_of_element(size_type i) {
505  return simplex_quality(gmm::index_ref_iterator
506  (pts.begin(), gmm::mat_col(t,i).begin()));
507  }
508 
509  void control_mesh_surface(void) {
510  mesh m;
511  adapt_mesh(m, 1);
512  dal::bit_vector ii = m.convex_index(), points_to_project;
513  size_type ic, ipt;
514  for (ic << ii; ic != size_type(-1); ic << ii) {
515  for (short_type f = 0; f <= N; ++f) {
516  if (!m.is_convex_having_neighbor(ic,f)) {
517  for (unsigned i = 0; i < N; ++i) {
518  ipt = m.ind_points_of_face_of_convex(ic, f)[i];
519  if (pts_attr[ipt]->constraints.card() == 0)
520  points_to_project.add(ipt);
521  else if ((*dist)(pts[ipt]) < -1e-2)
522  cout << "WARNING, point " << ipt
523  << " incoherent !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
524  }
525  }
526  }
527  }
528  if (points_to_project.card()) {
529  iter_wtcc = 0;
530  if (noisy > 1)
531  cout << "points to project : " << points_to_project << endl;
532  ii = points_to_project;
533  for (ipt << ii; ipt != size_type(-1); ipt << ii)
534  surface_projection_and_update_constraints(ipt);
535  }
536  }
537 
538  void suppress_flat_boundary_elements(void) {
539  size_type nb_deleted = 0;
540  do {
541  mesh m;
542  adapt_mesh(m, 1);
543  std::vector<base_small_vector> normals;
544  dal::bit_vector ii = m.convex_index();
545  dal::bit_vector flat_bound_elemnts;
546  size_type ic;
547 
548  for (ic << ii; ic != size_type(-1); ic << ii) {
549  scalar_type max_flatness = -2.0;
550  normals.resize(0);
551  for (short_type f = 0; f <= N; ++f) {
552  if (!m.is_convex_having_neighbor(ic,f)) {
553  if (quality_of_element(ic) < 1E-8) max_flatness = 1E-8;
554  else {
555  base_small_vector n = m.normal_of_face_of_convex(ic, f);
556  normals.push_back(n / gmm::vect_norm2(n));
557  }
558  }
559  }
560 
561  if (noisy > 1 && max_flatness != -2.0)
562  cout << "flatness of element " << ic << " : " << max_flatness;
563  if (normals.size() >= 2) {
564  if (noisy > 1) cout << "flatness of element " << ic << " : ";
565 
566  for (unsigned i = 1; i < normals.size(); ++i)
567  for (unsigned j = 0; j < i; ++j) {
568  scalar_type flatness=1.0-gmm::vect_sp(normals[i], normals[j]);
569  max_flatness = std::max(max_flatness, flatness);
570  if (noisy > 1) cout << flatness << " ";
571  }
572  }
573  if (max_flatness < boundary_threshold_flatness
574  && max_flatness!=-2.0) {
575  flat_bound_elemnts.add(ic);
576  if (noisy > 1) cout << " -> deleting";
577  }
578  if ((normals.size() >= 2 || max_flatness != -2.0) && noisy > 1)
579  cout << endl;
580 
581  }
582  ii = flat_bound_elemnts; nb_deleted = flat_bound_elemnts.card();
583  for (ic = ii.take_last(); ic != size_type(-1); ic = ii.take_last())
584  delete_element(ic);
585  } while (nb_deleted > 0);
586  }
587 
588 
589  bool try_projection(base_node &X)
590  { return getfem::try_projection(*dist, X); }
591 
592  void projection(base_node &X) {
593  base_small_vector G(X.size());
594  scalar_type d = dist->grad(X, G);
595  size_type it(0);
596  if (d > 0.0)
597  while (gmm::abs(d) > 1e-10) {
598  ++it;
599  GMM_ASSERT1(it <= 10000, "Object empty, or bad signed distance");
600 // cout << "iter " << it << " X = " << X << " dist = " << d <<
601 // " grad = " << G << endl;
602  gmm::add(gmm::scaled(G, -d / gmm::vect_norm2_sqr(G)), X);
603  d = dist->grad(X, G);
604  }
605  }
606 
607  void surface_projection(base_node &X) {
608  base_small_vector G;
609  scalar_type d = dist->grad(X, G);
610  size_type it(0);
611  while (gmm::abs(d) > 1e-10) {
612  ++it;
613  GMM_ASSERT1(it <= 10000, "Object empty, or bad signed distance X="
614  << X << ", G=" << G << " d = " << d);
615 // if (it > 9980)
616 // cout << "iter " << it << " X = " << X << " dist = " << d <<
617 // " grad = " << G << endl;
618  gmm::add(gmm::scaled(G, -d / gmm::vect_norm2_sqr(G)), X);
619  d = dist->grad(X, G);
620  }
621  }
622 
623  void projection(base_node &X, dal::bit_vector& bv)
624  { projection(X); bv.clear(); (*dist)(X,bv); }
625 
626  void constraint_projection(base_node &X, size_type cnum) {
627  base_small_vector G;
628  scalar_type d = constraints[cnum]->grad(X, G);
629  while (gmm::abs(d) > 1e-10) {
630  gmm::add(gmm::scaled(G, -d / gmm::vect_norm2_sqr(G)), X);
631  d=constraints[cnum]->grad(X, G);
632  }
633  }
634 
635  bool pure_multi_constraint_projection(base_node &X,
636  const dal::bit_vector &cts) {
637  getfem::pure_multi_constraint_projection(constraints, X, cts);
638 
639 // base_node oldX;
640 // size_type cnt = 0;
641 // do {
642 // oldX = X;
643 // for (dal::bv_visitor ic(cts); !ic.finished(); ++ic)
644 // constraint_projection(X,ic);
645 // ++cnt;
646 // } while (cts.card() && gmm::vect_dist2(oldX,X) > 1e-14 && cnt < 1000);
647 // if (cnt >= 1000) return false;
648 // if (cts.card()) {
649  dal::bit_vector ct2;
650  (*dist)(X, ct2);
651  return ct2.contains(cts);
652 // }
653 // return true;
654  }
655 
656  bool multi_constraint_projection(base_node &X,
657  const dal::bit_vector &cts) {
658  if (!(cts.card())) { projection(X); return true; }
659  else {
660  base_node oldX;
661  size_type cnt = 0;
662  do {
663  oldX = X;
664  for (dal::bv_visitor ic(cts); !ic.finished(); ++ic)
665  constraint_projection(X,ic);
666  projection(X);
667  ++cnt;
668  } while (gmm::vect_dist2(oldX,X) > 1e-14 && cnt < 1000);
669  if (cnt >= 1000) return false;
670  dal::bit_vector ct2;
671  (*dist)(X, ct2);
672  return ct2.contains(cts);
673  }
674  }
675 
676  void tangential_displacement(base_small_vector &X,
677  const dal::bit_vector &cts) {
678  base_small_vector G;
679  base_matrix normals(N, cts.card());
680  size_type cnt = 0;
681  for (dal::bv_visitor ic(cts); !ic.finished(); ++ic) {
682  constraints[ic]->grad(X, G);
683  gmm::copy(G, gmm::mat_col(normals, cnt));
684  for (size_type k=0; k < cnt; ++k) {
685  gmm::add(gmm::scaled(gmm::mat_col(normals, k),
686  -gmm::vect_sp(gmm::mat_col(normals,k),
687  gmm::mat_col(normals,cnt))),
688  gmm::mat_col(normals,cnt));
689  scalar_type n = gmm::vect_norm2(gmm::mat_col(normals, cnt));
690  if (n < 1e-8) continue;
691  gmm::scale(gmm::mat_col(normals, cnt), 1./n);
692  gmm::add(gmm::scaled(gmm::mat_col(normals, cnt),
693  -gmm::vect_sp(X,gmm::mat_col(normals, cnt))),X);
694  ++cnt;
695  }
696  }
697  }
698 
699  const pt_attribute* get_attr(bool fixed, const dal::bit_vector &bv) {
700  pt_attribute a;
701  a.fixed = fixed;
702  a.constraints = bv;
703  return &(*attributes_set.insert(a).first);
704  }
705 
706  void surface_projection_and_update_constraints(size_type ip) {
707  surface_projection(pts[ip]);
708  dal::bit_vector new_cts;
709  (*dist)(pts[ip], new_cts);
710  pts_attr[ip] = get_attr(pts_attr[ip]->fixed, new_cts);
711  }
712 
713  void project_and_update_constraints(size_type ip) {
714  const dal::bit_vector& cts = pts_attr[ip]->constraints;
715  dal::bit_vector new_cts;
716  multi_constraint_projection(pts[ip], cts);
717  (*dist)(pts[ip], new_cts);
718  if (noisy > 1 && !new_cts.contains(cts)) {
719  cout << "Point #" << ip << " has been downgraded from "
720  << cts << " to " << new_cts << endl;
721  }
722  else if (noisy > 1 && new_cts.card() > cts.card()) {
723  cout << "Point #" << ip << " has been upgraded from "
724  << cts << " to " << new_cts << endl;
725  }
726  if (new_cts != cts) {
727  pts_attr[ip] = get_attr(pts_attr[ip]->fixed, new_cts);
728  iter_wtcc = 0;
729  }
730 
731  }
732 
733  template <class VECT> void move_carefully(size_type ip, const VECT &VV) {
734  base_node V(N); gmm::copy(VV, V);
735 // if (pts_attr[ip]->constraints.card() != 0) {
736 // base_small_vector grad;
737 // dist->grad(pts[ip], grad);
738 // scalar_type r = gmm::vect_norm2_sqr(grad);
739 // gmm::add(gmm::scaled(grad, -gmm::vect_sp(V, grad) / r), V);
740 // }
741  scalar_type norm = gmm::vect_norm2(V);
742  if (norm <= h0/scalar_type(4))
743  gmm::add(V, pts[ip]);
744  else
745  gmm::add(gmm::scaled(V, h0 / (scalar_type(4) * norm)), pts[ip]);
746  project_and_update_constraints(ip);
747  }
748 
749  template <class VECT> void move_carefully(const VECT &V) {
750  scalar_type norm_max(0), lambda(1);
751  size_type npt = gmm::vect_size(V) / N;
752  for (size_type i = 0; i < npt; ++i)
753  norm_max = std::max(norm_max, gmm::vect_norm2
754  (gmm::sub_vector(V, gmm::sub_interval(i*N, N))));
755  if (norm_max > h0/scalar_type(3.7))
756  lambda = h0 / (scalar_type(3.7) * norm_max);
757 
758  for (size_type i = 0; i < npt; ++i)
759  move_carefully(i, gmm::scaled(gmm::sub_vector
760  (V, gmm::sub_interval(i*N, N)),lambda));
761  }
762 
763  void distribute_points_regularly(const std::vector<base_node>
764  &fixed_points) {
765  size_type nbpt = 1;
766  std::vector<size_type> gridnx(N);
767  mesh m;
768  base_node eff_boxmin(N), eff_boxmax(N);
769  bool eff_box_init = false;
770 
771  for (size_type i=0; i < N; ++i)
772  h0 = std::min(h0, bounding_box_max[i] - bounding_box_min[i]);
773  GMM_ASSERT1(h0 >= 1E-10, "h0 = " << h0 << " too small, aborting.");
774 
775  for (size_type i=0; i < N; ++i) {
776  scalar_type h = h0;
777  if (N == 2 && i == 1) h = sqrt(3.)/2. * h0;
778  gridnx[i]=1+(size_type)((bounding_box_max[i]-bounding_box_min[i])/h);
779  nbpt *= gridnx[i];
780  }
781 
782  /* build the regular grid and filter points outside */
783  for (size_type i=0; i < fixed_points.size(); ++i) {
784  if ((*dist)(fixed_points[i]) < geps &&
785  m.search_point(fixed_points[i]) == size_type(-1)) {
786  m.add_point(fixed_points[i]);
787  pts.push_back(fixed_points[i]);
788  pts_attr.push_back(get_attr(true,dal::bit_vector()));
789  }
790  else if (noisy > 0)
791  cout << "Removed duplicate fixed point: "<<fixed_points[i]<<"\n";
792  }
793  base_node P(N), Q(N);
794  for (size_type i=0; i < nbpt; ++i) {
795  for (size_type k=0, r = i; k < N; ++k) {
796  unsigned p = unsigned(r % gridnx[k]);
797  P[k] = p * (bounding_box_max[k] - bounding_box_min[k]) /
798  scalar_type((gridnx[k]-1)) + bounding_box_min[k];
799  if (N==2 && k==0 && ((r/gridnx[0])&1)==1) P[k] += h0/2;
800  r /= gridnx[k];
801  }
802 
803  dal::bit_vector co;
804  if ((prefind == 1 && (*dist)(P) < 0) || prefind == 2) {
805  for (size_type k = 0; k < constraints.size() && co.card() < N; ++k) {
806  gmm::copy(P, Q);
807  if (gmm::abs((*(constraints[k]))(Q)) < h0) {
808  constraint_projection(Q, k);
809  if ((*dist)(Q) > -geps &&
810  (gmm::vect_dist2(P, Q) < h0 / scalar_type(2))) co.add(k);
811  }
812  }
813  }
814  gmm::copy(P, Q);
815  if (co.card() > 0) {
816  bool ok = pure_multi_constraint_projection(Q, co);
817  if (!ok || gmm::abs((*dist)(Q)) > geps) { gmm::copy(P, Q); co.clear(); }
818  }
819 
820  if (prefind == 3) {
821  try_projection(Q);
822  }
823 
824  if ((*dist)(Q) < geps) {
825  if (m.search_point(Q) == size_type(-1)) {
826  //cout << "adding point : " << Q << endl;
827  if (!eff_box_init)
828  { eff_boxmin = eff_boxmax = Q; eff_box_init = true; }
829  else for (size_type k = 0; k < N; ++k) {
830  eff_boxmin[k] = std::min(eff_boxmin[k], Q[k]);
831  eff_boxmax[k] = std::max(eff_boxmax[k], Q[k]);
832  }
833  m.add_point(Q); pts.push_back(Q);
834  pts_attr.push_back(get_attr(false, co));
835  }
836  }
837  }
838  if (noisy > 0) cout << "effective bounding box : " << eff_boxmin << " : " << eff_boxmax << endl;
839  if (prefind == 3) {
840  h0 = std::min(h0, gmm::vect_dist2(eff_boxmin, eff_boxmax)
841  / scalar_type(2));
842  }
843  }
844 
845  void add_point_hull(void) {
846  if (dist_point_hull > 0) {
847  size_type nbpt = pts.size(), nbadd(0);
848  base_node P, Q;
849  base_small_vector V;
850  for (unsigned i=0; i < nbpt; ++i) {
851  if (pts_attr[i]->constraints.card()) {
852  P = pts[i];
853  dist->grad(P, V);
854  scalar_type d = gmm::vect_norm2(V);
855  if (d > 0) {
856  P += V * (dist_point_hull*h0/d);
857  if ((*dist)(P)*sqrt(scalar_type(N)) > dist_point_hull*h0) {
858  Q = P;
859  projection(Q);
860  if (gmm::vect_dist2(P, Q) > dist_point_hull*h0/scalar_type(2))
861  { pts.push_back(P); ++nbadd; }
862  }
863  }
864  }
865  }
866  if (noisy > 1) cout << "point hull: " << nbadd << " points added\n";
867  }
868  }
869 
870 
871  scalar_type pts_dist_max(const std::vector<base_node> &A,
872  const std::vector<base_node> &B) {
873  scalar_type dist_max = 0;
874  for (size_type i=0; i < pts.size(); ++i)
875  dist_max = std::max(dist_max, gmm::vect_dist2_sqr(A[i],B[i]));
876  return sqrt(dist_max);
877  }
878 
879  struct cleanup_points_compare {
880  const std::vector<base_node> &pts;
881  const std::vector<const pt_attribute*> &attr;
882  cleanup_points_compare(const std::vector<base_node> &pts_,
883  const std::vector<const pt_attribute*> &attr_)
884  : pts(pts_), attr(attr_) {}
885  bool operator()(size_type a, size_type b) {
886  if (attr[a] != attr[b]) return attr[a] < attr[b];
887  return pts[a] < pts[b];
888  }
889  };
890  void cleanup_points() {
891  std::vector<size_type> idx(pts.size());
892  for (size_type i=0; i < idx.size(); ++i) idx[i] = i;
893 
894  std::sort(idx.begin(), idx.end(), cleanup_points_compare(pts,pts_attr));
895  bgeot::kdtree tree;
896  bgeot::kdtree_tab_type neighbors;
897  dal::bit_vector keep_pts; keep_pts.add(0,idx.size());
898  for (size_type i=0, i0=0; i < idx.size(); ++i) {
899  const base_node &P = pts[idx[i]];
900  const pt_attribute *a = pts_attr[idx[i]];
901  tree.add_point_with_id(P,i);
902  if (i == idx.size()-1 || a != pts_attr[idx[i+1]]) {
903  for (size_type j=i0; j < i+1; ++j) {
904  base_node bmin = P, bmax = P;
905  scalar_type h = h0*edge_len(P);
906  for (size_type k = 0; k < N; ++k)
907  { bmin[k] -= h/20.; bmax[k] += h/20.; }
908 
909  tree.points_in_box(neighbors, bmin, bmax);
910  for (size_type k=0; k < neighbors.size(); ++k) {
911  if (neighbors[k].i != i && keep_pts.is_in(neighbors[k].i)
912  && keep_pts.is_in(i)) {
913  if (noisy > 0)
914  cout << "point #" << i << " " << P
915  << " is too near from point #" << neighbors[k].i
916  << pts[idx[neighbors[k].i]] << " : will be removed\n";
917  keep_pts.sup(i);
918  }
919  }
920  }
921  tree.clear();
922  i0 = i+1;
923  }
924  }
925  pts_prev.resize(keep_pts.card());
926  size_type cnt = 0;
927  std::vector<const pt_attribute*> pts_attr2(keep_pts.card());
928  for (dal::bv_visitor i(keep_pts); !i.finished(); ++i, ++cnt) {
929  pts_prev[cnt].swap(pts[idx[i]]);
930  pts_attr2[cnt] = pts_attr[idx[i]];
931  }
932  pts_attr.swap(pts_attr2);
933  pts.resize(pts_prev.size());
934  std::copy(pts_prev.begin(), pts_prev.end(), pts.begin());
935  }
936 
937  scalar_type worst_q;
938  base_node worst_q_P;
939 
940  void select_elements(int version) {
941  size_type nbpt = pts.size();
942 
943  worst_q = 1.;
944  base_node weights(N+1);
945  for (size_type i=0; i < gmm::mat_ncols(t); ) {
946  bool ext_simplex = false;
947  // bool boundary_simplex = true;
948  bool on_boundary_simplex = false;
949  bool is_bridge_simplex = false;
950  scalar_type q(0), dG(0);
951  base_node G;
952 
953  for (size_type k=0; k <= N; ++k)
954  if (t(k, i) >= nbpt) ext_simplex = true;
955 
956  if (!ext_simplex) {
957  G = pts[t(0,i)];
958  for (size_type k=1; k <= N; ++k) G += pts[t(k,i)];
959  gmm::scale(G, scalar_type(1)/scalar_type(N+1));
960  dG = (*dist)(G);
961  gmm::clear(weights);
962 
963  q = quality_of_element(i);
964 
965  for (size_type k=0; k <= N; ++k) {
966  if (!(pts_attr[t(k,i)]->constraints.card() == 0))
967  on_boundary_simplex = true;
968  // else
969  // boundary_simplex = false;
970  }
971 
972  if (version == 1 && on_boundary_simplex)
973  for (size_type k=1; k < N+1; ++k)
974  for (size_type l=0; l < k; ++l) {
975  dal::bit_vector all_cts = pts_attr[t(k,i)]->constraints
976  | pts_attr[t(l,i)]->constraints;
977  if (/* gmm::vect_dist2(pts[t(k,i)], pts[t(l,i)]) > h0 && */
978  !(pts_attr[t(k,i)]->constraints.contains(all_cts))
979  && !(pts_attr[t(l,i)]->constraints.contains(all_cts))
980  && (*dist)(0.5*(pts[t(k,i)] + pts[t(l,i)])) > 0.)
981  is_bridge_simplex = true;
982  }
983  }
984  if (ext_simplex || dG > 0 || is_bridge_simplex || q < 1e-14) {
985  delete_element(i);
986  } else {
987  ++i;
988  if (q < worst_q) {
989  worst_q = q;
990  worst_q_P = G*(scalar_type(1)/scalar_type(N+1));
991  }
992  }
993  }
994 
995  }
996 
997  void adapt_mesh(mesh &m, size_type degree) {
998  std::vector<base_node> cvpts(N+1), cvpts2;
999  size_type cvnum;
1000  m.clear();
1001  for (size_type ip=0; ip < pts.size(); ++ip) {
1002  size_type z;
1003  bgeot::small_vector<scalar_type> P = pts[ip];
1004  while ((z = m.add_point(P)) != ip) {
1005  if (noisy > 0) cout << "WARNING : points are too near ...\n";
1006  bgeot::small_vector<scalar_type> Z(N); gmm::fill_random(Z);
1007  gmm::add(gmm::scaled(Z, h0/1000.0), P);
1008  }
1009  }
1010  for (size_type i=0; i < t.size()/(N+1); ++i) {
1011  for (size_type k=0; k < N+1; ++k) cvpts[k] = pts[t[i*(N+1)+k]];
1012  if (degree == 1) {
1013  cvnum = m.add_convex(bgeot::simplex_geotrans(N,1), &t[i*(N+1)]);
1014  assert(cvnum == i);
1015  } else {
1016  bgeot::pgeometric_trans pgt =
1017  bgeot::simplex_geotrans(N,short_type(degree));
1018  cvpts2.resize(pgt->nb_points());
1019  for (size_type k=0; k < pgt->nb_points(); ++k) {
1020  cvpts2[k] = bgeot::simplex_geotrans(N,1)->transform
1021  (pgt->convex_ref()->points()[k], cvpts);
1022  }
1023  cvnum = m.add_convex_by_points(pgt, cvpts2.begin());
1024  assert(cvnum == i);
1025  }
1026  }
1027  if (degree>1) {
1028  //m.optimize_structure();
1029  getfem::mesh_region border_faces;
1030  getfem::outer_faces_of_mesh(m, border_faces);
1031  dal::bit_vector ptdone; ptdone.sup(0,m.points_index().last_true());
1032  for (getfem::mr_visitor it(border_faces); !it.finished(); ++it) {
1033  mesh::ind_pt_face_ct fpts_
1034  = m.ind_points_of_face_of_convex(it.cv(), it.f());
1035  std::vector<size_type> fpts(fpts_.size());
1036  std::copy(fpts_.begin(), fpts_.end(), fpts.begin());
1037  interpolate_face(m, ptdone, fpts,
1038  m.trans_of_convex(it.cv())->structure()
1039  ->faces_structure()[it.f()]);
1040  }
1041  }
1042  }
1043 
1044 
1045 
1046  void interpolate_face(mesh &m, dal::bit_vector& ptdone,
1047  const std::vector<size_type>& ipts,
1048  bgeot::pconvex_structure cvs) {
1049  if (cvs->dim() == 0) return;
1050  else if (cvs->dim() > 1) {
1051  std::vector<size_type> fpts;
1052  for (short_type f=0; f < cvs->nb_faces(); ++f) {
1053  fpts.resize(cvs->nb_points_of_face(f));
1054  for (size_type k=0; k < fpts.size(); ++k)
1055  fpts[k] = ipts[cvs->ind_points_of_face(f)[k]];
1056  interpolate_face(m,ptdone,fpts,cvs->faces_structure()[f]);
1057  }
1058  }
1059  dal::bit_vector cts; size_type cnt = 0;
1060  for (size_type i=0; i < ipts.size(); ++i) {
1061  if (ipts[i] < pts.size()) {
1062  if (cnt == 0) cts = pts_attr[ipts[i]]->constraints;
1063  else cts &= pts_attr[ipts[i]]->constraints;
1064  ++cnt;
1065  }
1066  }
1067  if (cts.card()) {
1068  // dal::bit_vector new_cts;
1069  for (size_type i=0; i < ipts.size(); ++i) {
1070  if (ipts[i] >= pts.size() && !ptdone[ipts[i]]) {
1071  base_node &P = m.points()[ipts[i]];
1072  multi_constraint_projection(P, cts);
1073  // (*dist)(P, new_cts);
1074  }
1075  }
1076  }
1077  }
1078 
1079  void special_constraints_management(void) {
1080 
1081  bgeot::kdtree tree;
1082  bgeot::kdtree_tab_type neighbors;
1083  bool tree_empty = true;
1084  mesh::ind_set iAneighbors, iBneighbors, common_pts;
1085 
1086  attractor_points.resize(0); attracted_points.resize(0);
1087 
1088  for (dal::bv_visitor ie(edges_mesh.convex_index());
1089  !ie.finished(); ++ie) {
1090  size_type iA = edges_mesh.ind_points_of_convex(ie)[0];
1091  size_type iB = edges_mesh.ind_points_of_convex(ie)[1];
1092  // if (L[ie] > L0[ie]) continue;
1093  if (pts_attr[iA] == pts_attr[iB] ||
1094  pts_attr[iA]->constraints.card() == 0 ||
1095  pts_attr[iB]->constraints.card() == 0) continue;
1096  if (pts_attr[iA]->constraints == pts_attr[iB]->constraints)
1097  continue;
1098  dal::bit_vector bv1(pts_attr[iA]->constraints);
1099  bv1.setminus(pts_attr[iB]->constraints);
1100  dal::bit_vector bv2(pts_attr[iB]->constraints);
1101  bv2.setminus(pts_attr[iA]->constraints);
1102  if (bv1.card() && bv2.card()) {
1103  bv1 |= bv2;
1104  edges_mesh.ind_points_to_point(iA, iAneighbors);
1105  edges_mesh.ind_points_to_point(iB, iBneighbors);
1106  common_pts.resize(0);
1107  for (size_type i = 0; i < iAneighbors.size(); ++i)
1108  if (std::find(iBneighbors.begin(), iBneighbors.end(),
1109  iAneighbors[i]) != iBneighbors.end())
1110  common_pts.push_back(iAneighbors[i]);
1111  bool do_projection = true;
1112  if ((*dist)(.5*(pts[iA]+pts[iB])) < 0) {
1113  for (mesh::ind_set::iterator it = common_pts.begin();
1114  it != common_pts.end(); ++it) {
1115  if (pts_attr[*it]->constraints.contains(bv1)) {
1116  do_projection = false;
1117  break;
1118  }
1119  }
1120  }
1121  if (do_projection) {
1122 
1123  if (pts_attr[iA]->constraints.card()
1124  < pts_attr[iB]->constraints.card()) std::swap(iA,iB);
1125 
1126  base_node PA = pts[iA];
1127  bool okA = multi_constraint_projection(PA, bv1);
1128  base_node PB = pts[iB];
1129  bool okB = multi_constraint_projection(PB, bv1);
1130 
1131  if (okB && !okA)
1132  { std::swap(PA,PB); std::swap(iA,iB); std::swap(okB, okA); }
1133 
1134  if (okB && (gmm::vect_dist2(PA,pts[iA])
1135  > 1.1*gmm::vect_dist2(PB,pts[iB]))) {
1136  // 1.5 au lieu de 1.1 ?
1137  std::swap(iA,iB); std::swap(PA,PB);
1138  }
1139 
1140 
1141  if (okA && gmm::vect_dist2(PA, pts[iA]) < h0*0.75) {
1142 
1143  if (tree_empty) {
1144  for (size_type i=0; i < pts.size(); ++i)
1145  tree.add_point_with_id(pts[i],i);
1146  tree_empty = false;
1147  }
1148 
1149  base_node bmin = PA, bmax = PA;
1150  for (size_type k = 0; k < N; ++k)
1151  { bmin[k] -= h0/1.8; bmax[k] += h0/1.8; }
1152  tree.points_in_box(neighbors, bmin, bmax);
1153  for (size_type k=0; k < neighbors.size(); ++k) {
1154  if (neighbors[k].i != iA && neighbors[k].i != iB)
1155  do_projection = false;
1156  }
1157 
1158  if (do_projection) {
1159  attractor_points.push_back(PA);
1160  attracted_points.push_back(iA);
1161  }
1162  }
1163  }
1164  }
1165  }
1166  }
1167 
1168 
1169  void running_delaunay(bool mct) {
1170  if (noisy > 0)
1171  cout << "NEW DELAUNAY, running on " << pts.size() << " points\n";
1172  size_type nbpt = pts.size();
1173  add_point_hull();
1174  bgeot::qhull_delaunay(pts, t);
1175  pts.resize(nbpt);
1176  if (noisy > 1) cout << "number of elements before selection = "
1177  << gmm::mat_ncols(t) << "\n";
1178  if (mct) {
1179  select_elements(0);
1180  edges_mesh.clear();
1181  for (size_type i=0; i < gmm::mat_ncols(t); ++i)
1182  for (size_type j=0; j < N+1; ++j)
1183  for (size_type k=j+1; k < N+1; ++k)
1184  edges_mesh.add_segment(t(j,i), t(k,i));
1185  special_constraints_management();
1186  }
1187  select_elements(1);
1188  if (noisy > 0) cout << "number of elements after selection = "
1189  << gmm::mat_ncols(t) << "\n";
1190  edges_mesh.clear();
1191  for (size_type i=0; i < gmm::mat_ncols(t); ++i)
1192  for (size_type j=0; j < N+1; ++j)
1193  for (size_type k=j+1; k < N+1; ++k)
1194  edges_mesh.add_segment(t(j,i), t(k,i));
1195  }
1196 
1197  void standard_move_strategy(base_vector &X) {
1198  for (dal::bv_visitor ie(edges_mesh.convex_index());
1199  !ie.finished(); ++ie) {
1200  size_type iA = edges_mesh.ind_points_of_convex(ie)[0];
1201  size_type iB = edges_mesh.ind_points_of_convex(ie)[1];
1202  base_node bar = pts[iB]-pts[iA];
1203  scalar_type F = std::max(L0[ie]-L[ie], 0.);
1204 
1205  if (F) {
1206  base_node Fbar = (bar)*(F/L[ie]);
1207 
1208  if (!pts_attr[iA]->fixed) // pts[iA] -= deltat*Fbar;
1209  gmm::add(gmm::scaled(Fbar, -deltat),
1210  gmm::sub_vector(X, gmm::sub_interval(iA*N, N)));
1211  if (!pts_attr[iB]->fixed) // pts[iB] += deltat*Fbar;
1212  gmm::add(gmm::scaled(Fbar, deltat),
1213  gmm::sub_vector(X, gmm::sub_interval(iB*N, N)));
1214  }
1215  }
1216  }
1217 
1218  void do_build_mesh(mesh &m,
1219  const std::vector<base_node> &fixed_points) {
1220 
1221  distribute_points_regularly(fixed_points);
1222  std::vector<base_node> pts2(pts.size(),base_node(N));
1223  size_type count = 0, count_id = 0, count_ct = 0;
1224  bool pt_changed = false;
1225  iter_wtcc = 0;
1226  attracted_points.resize(0);
1227  attractor_points.resize(0);
1228 
1229  do {
1230  if (noisy > 1) {
1231  cout << "Iter " << count << " / " << count + iter_max - iter_wtcc;
1232  if (count && pts_prev.size() == pts.size())
1233  cout << ", dist_max since last delaunay ="
1234  << pts_dist_max(pts, pts_prev) << ", tol=" << ttol*h0;
1235  cout << endl;
1236  }
1237  if (count==0 || pts_prev.size() != pts.size()
1238  || (pts_dist_max(pts, pts_prev) > ttol*h0 && count_id >= 5)) {
1239  size_type nbpt = pts.size();
1240  cleanup_points(); /* and copy pts to pts_prev */
1241  if (noisy == 1) cout << "Iter " << count << " ";
1242  bool mct = false;
1243  if (count_ct >= 20 || pt_changed) {
1244  mct = true;
1245  count_ct = 0;
1246  }
1247  running_delaunay(mct);
1248  pt_changed = nbpt != pts.size();
1249  count_id = 0;
1250  }
1251  ++count_id; ++count_ct;
1252  pts2 = pts;
1253 
1254  // computation of L and L0.
1255  size_type nbcv = edges_mesh.convex_index().card();
1256  GMM_ASSERT1(nbcv != 0, "no more edges!");
1257  L.resize(nbcv); L0.resize(nbcv);
1258  scalar_type sL = 0, sL0 = 0;
1259  for (dal::bv_visitor ie(edges_mesh.convex_index());
1260  !ie.finished(); ++ie) {
1261  const base_node &A = pts[edges_mesh.ind_points_of_convex(ie)[0]];
1262  const base_node &B = pts[edges_mesh.ind_points_of_convex(ie)[1]];
1263  base_node C(A); C+=B; C /= scalar_type(2);
1264  L[ie] = gmm::vect_dist2(A, B);
1265  L0[ie] = edge_len(C);
1266  sL += pow(L[ie],scalar_type(N));
1267  sL0 += pow(L0[ie],scalar_type(N));
1268  }
1269  gmm::scale(L0, L0mult * pow(sL/sL0, scalar_type(1)/scalar_type(N)));
1270 
1271  // Moving the points with standard strategy
1272  base_vector X(pts.size() * N);
1273  standard_move_strategy(X);
1274  for (size_type i = 0; i < attracted_points.size(); ++i) {
1275  size_type npt = attracted_points[i];
1276  gmm::copy(attractor_points[i] - pts[npt],
1277  gmm::sub_vector(X, gmm::sub_interval(npt*N, N)));
1278  }
1279 
1280  move_carefully(X);
1281 
1282  scalar_type maxdp = pts_dist_max(pts,pts2);
1283  if (noisy > 1)
1284  cout << ", maxdp = " << maxdp << ", ptol = "
1285  << ptol << " CV=" << sqrt(maxdp)*deltat/h0 << "\n";
1286  ++count; ++iter_wtcc;
1287 
1288  if (iter_wtcc == 100) control_mesh_surface();
1289 
1290  // m.clear();
1291  // for (size_type i=0; i < t.size()/(N+1); ++i)
1292  // m.add_convex_by_points(bgeot::simplex_geotrans(N,1),
1293  // dal::index_ref_iterator(pts.begin(), &t[i*(N+1)]));
1294  // char s[50]; sprintf(s, "toto%02d.mesh", count);
1295  // m.write_to_file(s);
1296 
1297  if ( (count > 40 && sqrt(maxdp)*deltat < ptol * h0)
1298  || iter_wtcc>iter_max || count > 10000) {
1299 
1300 // {
1301 // m.clear();
1302 // adapt_mesh(m,K);
1303 // m.optimize_structure();
1304 // getfem::vtk_export exp("toto1.vtk");
1305 // exp.exporting(m);
1306 // exp.write_mesh_quality(m);
1307 // }
1308 
1309  control_mesh_surface();
1310  size_type nbpt = pts.size();
1311  add_point_hull();
1312  bgeot::qhull_delaunay(pts, t);
1313  pts.resize(nbpt);
1314  select_elements((prefind == 3) ? 0 : 1);
1315  suppress_flat_boundary_elements();
1316 
1317 // {
1318 // m.clear();
1319 // adapt_mesh(m,K);
1320 // m.optimize_structure();
1321 // getfem::vtk_export exp("toto2.vtk");
1322 // exp.exporting(m);
1323 // exp.write_mesh_quality(m);
1324 // }
1325 
1326  if (prefind != 3) optimize_quality();
1327 
1328  // ajout d'un point au barycentre des elements trop plats :
1329 
1330  if (prefind != 3)
1331  for (unsigned cv = 0; cv < gmm::mat_ncols(t); ++cv) {
1332 
1333  if (quality_of_element(cv) < 0.05) {
1334  base_node G = pts[t(0,cv)];
1335  for (size_type k=1; k <= N; ++k) G += pts[t(k,cv)];
1336  gmm::scale(G, scalar_type(1)/scalar_type(N+1));
1337  pts.push_back(G);
1338  pts_attr.push_back(get_attr(false, dal::bit_vector()));
1339  }
1340  }
1341 
1342  if (pts.size() != nbpt) {
1343  control_mesh_surface();
1344  nbpt = pts.size();
1345  add_point_hull();
1346  bgeot::qhull_delaunay(pts, t);
1347  pts.resize(nbpt);
1348  select_elements((prefind == 3) ? 0 : 1);
1349  suppress_flat_boundary_elements();
1350 
1351 // {
1352 // m.clear();
1353 // adapt_mesh(m,K);
1354 // m.optimize_structure();
1355 // getfem::vtk_export exp("toto3.vtk");
1356 // exp.exporting(m);
1357 // exp.write_mesh_quality(m);
1358 // }
1359 
1360  if (prefind != 3) optimize_quality();
1361  }
1362  break;
1363  }
1364 
1365 
1366  } while (true);
1367 
1368  m.clear();
1369  adapt_mesh(m,K);
1370  // m.write_to_file("toto.mesh");
1371  m.optimize_structure();
1372 
1373 
1374 // getfem::vtk_export exp("toto4.vtk");
1375 // exp.exporting(m);
1376 // exp.write_mesh_quality(m);
1377 
1378 // getfem::stored_mesh_slice sl;
1379 // sl.build(m, getfem::slicer_explode(0.8), 8);
1380 // getfem::vtk_export exp2("totoq.vtk");
1381 // exp2.exporting(sl);
1382 // exp2.write_mesh();
1383 // exp2.write_mesh_quality(m);
1384 // getfem::dx_export exp3("totoq.dx");
1385 // exp3.exporting(sl);
1386 // exp3.write_mesh();
1387 
1388 
1389 
1390 
1391 // getfem::stored_mesh_slice slb; slb.build(m, getfem::slicer_boundary(m), 4);
1392 // getfem::stored_mesh_slice sl2;
1393 // getfem::mesh_slicer ms(m);
1394 
1395 
1396 // getfem::slicer_build_stored_mesh_slice bb(sl2);
1397 // ms.push_back_action(bb);
1398 // getfem::convex_face_ct cvlst;
1399 // for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
1400 // scalar_type q = m.convex_quality_estimate(cv);
1401 // if (q< 0.2)
1402 // cvlst.push_back(getfem::convex_face(cv));
1403 // //cout << "cv " << cv << ": q=" << q << "\n";
1404 // }
1405 // //ms.exec(3, cvlst);
1406 // sl2.merge(slb);
1407 
1408 
1409 // getfem::vtk_export exp3("totoq2.vtk");
1410 // exp3.exporting(sl2);
1411 // exp3.write_mesh();
1412 // exp3.write_mesh_quality(m);
1413  }
1414 
1415 
1416  };
1417 
1418 
1419 
1420  void build_mesh(mesh &m, const pmesher_signed_distance& dist,
1421  scalar_type h0, const std::vector<base_node> &fixed_points,
1422  size_type K, int noise, size_type iter_max, int prefind,
1423  scalar_type dist_point_hull,
1424  scalar_type boundary_threshold_flatness) {
1425  mesher mg(K, dist, getfem::mvf_constant(1), h0, m, fixed_points, noise,
1426  iter_max, prefind, dist_point_hull, boundary_threshold_flatness);
1427  }
1428 
1429 }
1430 
gmm::lu_factor
size_type lu_factor(DenseMatrix &A, lapack_ipvt &ipvt)
LU Factorization of a general (dense) matrix (real or complex).
Definition: gmm_dense_lu.h:129
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
bgeot::kdtree_tab_type
std::vector< index_node_pair > kdtree_tab_type
store a set of points with associated indexes.
Definition: bgeot_kdtree.h:68
gmm::is_hermitian
bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol=magnitude_of_linalg(MAT)(-1))
*‍/
Definition: gmm_blas.h:2229
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
getfem_mesher.h
An experimental mesher.
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:597
bgeot::mesh_structure::ind_points_of_convex
const ind_set & ind_points_of_convex(size_type ic) const
Return a container to the list of points attached to convex ic.
Definition: bgeot_mesh_structure.h:107
bgeot::mesh_structure::clear
void clear()
erase the mesh
Definition: bgeot_mesh_structure.cc:205
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
bgeot::mesh_structure::convex_index
const dal::bit_vector & convex_index() const
Return the list of valid convex IDs.
Definition: bgeot_mesh_structure.h:91
bgeot::kdtree
Balanced tree over a set of points.
Definition: bgeot_kdtree.h:102
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:237
bgeot::mesh_structure::ind_points_to_point
void ind_points_to_point(size_type, ind_set &) const
Return a container of the points attached (via an edge) to point ip.
Definition: bgeot_mesh_structure.cc:159
bgeot::kdtree::clear
void clear()
reset the tree, remove all points
Definition: bgeot_kdtree.h:113
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
bgeot::small_vector< scalar_type >
gmm::fill_random
void fill_random(L &l, double cfill)
*‍/
Definition: gmm_blas.h:154
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
bgeot::mesh_structure
Mesh structure definition.
Definition: bgeot_mesh_structure.h:73
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:54
gmm::mat_trace
linalg_traits< M >::value_type mat_trace(const M &m)
Trace of a matrix.
Definition: gmm_blas.h:528
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::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
getfem::outer_faces_of_mesh
void APIDECL outer_faces_of_mesh(const mesh &m, const dal::bit_vector &cvlst, convex_face_ct &flist)
returns a list of "exterior" faces of a mesh (i.e.
Definition: getfem_mesh.cc:822
bgeot::equilateral_simplex_of_reference
pconvex_ref equilateral_simplex_of_reference(dim_type nc)
equilateral simplex (degree 1).
Definition: bgeot_convex_ref.cc:831
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:120