GetFEM  5.4.2
bgeot_node_tab.cc
1 /*===========================================================================
2 
3  Copyright (C) 2007-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 
23 #include "getfem/bgeot_node_tab.h"
24 
25 namespace bgeot {
26 
27  bool node_tab::component_comp::operator()(size_type i1,size_type i2) const {
28  if (i1 == i2) return false;
29  const base_node &pt1((i1 == size_type(-1)) ? *c : (*vbn)[i1]);
30  const base_node &pt2((i2 == size_type(-1)) ? *c : (*vbn)[i2]);
31  unsigned d = pt1.size();
32  base_small_vector::const_iterator it = v.begin();
33  base_node::const_iterator it1 = pt1.begin(), it2 = pt2.begin();
34  scalar_type a(0);
35  for (size_type i = 0; i < d; ++i) a += (*it++) * (*it1++ - *it2++);
36  if (a != scalar_type(0)) return a < 0;
37  if (i1 == size_type(-1) || i2 == size_type(-1)) return false;
38  return i1 < i2;
39  }
40 
41  node_tab::component_comp::component_comp
42  (const dal::dynamic_tas<base_node> &vbn_, const base_node &c_, unsigned dim)
43  : vbn(&vbn_), c(&c_), v(dim) {
44  do gmm::fill_random(v); while (gmm::vect_norm2(v) == 0);
45  gmm::scale(v, scalar_type(1) / gmm::vect_norm2(v) );
46  }
47 
48  void node_tab::add_sorter(void) const {
49  if (sorters.size() > 1)
50  GMM_WARNING3("Multiple sort needed for node tab : " << sorters.size()+1);
51  sorters.push_back(sorter(component_comp(*this, c, dim_)));
52  for (dal::bv_visitor i(index()); !i.finished(); ++i)
53  sorters.back().insert(size_type(i));
54  }
55 
56  size_type node_tab::search_node(const base_node &pt,
57  const scalar_type radius) const {
58  if (card() == 0 || radius < 0.)
59  return size_type(-1);
60 
61  scalar_type eps_radius = std::max(eps, radius);
62  for (size_type is = 0; is < 5; ++is) {
63  if (is >= sorters.size()) add_sorter();
64 
65  c = pt - eps_radius * sorters[is].key_comp().v;
66 
67  sorter::const_iterator it = sorters[is].lower_bound(size_type(-1));
68  scalar_type up_bound
69  = gmm::vect_sp(pt, sorters[is].key_comp().v) + 2*eps_radius;
70  size_type count = 0;
71  // if (is > 0) cout << "begin loop " << " v = " << sorters[is].key_comp().v << "sp c = " << gmm::vect_sp(c, sorters[is].key_comp().v) << " eps_radius = " << eps_radius << " max_radius " << max_radius << endl;
72  for (; it != sorters[is].end() && count < 20; ++it, ++count) {
73 
74  const base_node &pt2 = (*this)[*it];
75 
76 // if (count > 0) {
77 // cout << "count " << count << " is = " << is << " pt = " << pt << " pt2 = " << pt2 << " sp = " << gmm::vect_sp(pt2, sorters[is].key_comp().v) << " spinit = " << gmm::vect_sp(pt, sorters[is].key_comp().v) << endl;
78 // }
79 
80  if (gmm::vect_dist2(pt, pt2) < eps_radius) return *it;
81  if (gmm::vect_sp(pt2, sorters[is].key_comp().v) > up_bound)
82  return size_type(-1);
83  }
84  if (it == sorters[is].end()) return size_type(-1);
85  }
86  GMM_ASSERT1(false, "Problem in node structure !!");
87  }
88 
89  void node_tab::clear() {
90  dal::dynamic_tas<base_node>::clear();
91  sorters = std::vector<sorter>();
92  max_radius = scalar_type(1e-60);
93  eps = max_radius * prec_factor;
94  }
95 
96  size_type node_tab::add_node(const base_node &pt,
97  const scalar_type radius,
98  bool remove_duplicated_nodes) {
99  scalar_type npt = gmm::vect_norm2(pt);
100  max_radius = std::max(max_radius, npt);
101  eps = max_radius * prec_factor;
102 
103  if (this->card() == 0)
104  dim_ = pt.size();
105  else
106  GMM_ASSERT1(dim_ == pt.size(), "Nodes should have the same dimension");
107  size_type id(-1);
108  if (remove_duplicated_nodes && radius >= 0.)
109  id = search_node(pt, radius);
110  if (id == size_type(-1)) {
111  id = dal::dynamic_tas<base_node>::add(pt);
112  for (size_type i = 0; i < sorters.size(); ++i) {
113  sorters[i].insert(id);
114  GMM_ASSERT3(sorters[i].size() == card(), "internal error");
115  }
116  }
117  return id;
118  }
119 
120  void node_tab::swap_points(size_type i, size_type j) {
121  if (i != j) {
122  bool existi = index().is_in(i), existj = index().is_in(j);
123  for (size_type is = 0; is < sorters.size(); ++is) {
124  if (existi) sorters[is].erase(i);
125  if (existj) sorters[is].erase(j);
126  }
127  dal::dynamic_tas<base_node>::swap(i, j);
128  for (size_type is = 0; is < sorters.size(); ++is) {
129  if (existi) sorters[is].insert(j);
130  if (existj) sorters[is].insert(i);
131  GMM_ASSERT3(sorters[is].size() == card(), "internal error");
132  }
133  }
134  }
135 
136  void node_tab::sup_node(size_type i) {
137  if (index().is_in(i)) {
138  for (size_type is = 0; is < sorters.size(); ++is) {
139  sorters[is].erase(i);
140  GMM_ASSERT3(sorters[is].size()+1 == card(), "Internal error");
141  // if (sorters[is].size()+1 != card()) { resort(); }
142  }
143  dal::dynamic_tas<base_node>::sup(i);
144 
145  }
146  }
147 
148  void node_tab::translation(const base_small_vector &V) {
149  for (dal::bv_visitor i(index()); !i.finished(); ++i) (*this)[i] += V;
150  resort();
151  }
152 
153  void node_tab::transformation(const base_matrix &M) {
154  base_small_vector w(M.nrows());
155  GMM_ASSERT1(gmm::mat_nrows(M) != 0 && gmm::mat_ncols(M) == dim(),
156  "invalid dimensions for the transformation matrix");
157  dim_ = unsigned(gmm::mat_nrows(M));
158  for (dal::bv_visitor i(index()); !i.finished(); ++i) {
159  w = (*this)[i];
160  gmm::resize((*this)[i], dim_);
161  gmm::mult(M,w,(*this)[i]);
162  }
163  resort();
164  }
165 
166  node_tab::node_tab(scalar_type prec_loose) {
167  max_radius = scalar_type(1e-60);
168  sorters.reserve(5);
169  prec_factor = gmm::default_tol(scalar_type()) * prec_loose;
170  eps = max_radius * prec_factor;
171  }
172 
173  node_tab::node_tab(const node_tab &t)
174  : dal::dynamic_tas<base_node>(t), sorters(), eps(t.eps),
175  prec_factor(t.prec_factor), max_radius(t.max_radius), dim_(t.dim_) {}
176 
177  node_tab &node_tab::operator =(const node_tab &t) {
178  dal::dynamic_tas<base_node>::operator =(t);
179  sorters = std::vector<sorter>();
180  eps = t.eps; prec_factor = t.prec_factor;
181  max_radius = t.max_radius; dim_ = t.dim_;
182  return *this;
183  }
184 
185 }
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
bgeot_node_tab.h
Structure which dynamically collects points identifying points that are nearer than a certain very sm...
dal
Dynamic Array Library.
Definition: dal_backtrace.cc:29
bgeot::small_vector< scalar_type >
gmm::fill_random
void fill_random(L &l, double cfill)
*‍/
Definition: gmm_blas.h:154
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:27