Delaunay.cpp
1 #include <fstream>
2 #include "Delaunay.hpp"
3 #include <vector>
4 #include <cmath>
5 #include "../misc/triplet.hpp"
6 #include <boost/foreach.hpp>
7 #include <boost/multiprecision/cpp_dec_float.hpp>
8 #ifdef RICH_MPI
9 #include <mpi.h>
10 #endif // RICH_MPI
11 
12 using namespace std;
13 
14 namespace
15 {
16 #ifdef RICH_MPI
17  void PeriodicGetRidDuplicatesSingle(std::vector<int> &duplicate, std::vector<Vector2D> &toadd)
18  {
19  if (duplicate.empty())
20  return;
21  assert(duplicate.size() == toadd.size());
22  std::vector<size_t> indeces = sort_index(duplicate);
23  toadd = VectorValues(toadd, indeces);
24  duplicate = VectorValues(duplicate, indeces);
25  size_t N = duplicate.size();
26  std::vector<double> dx(N, 0);
27  std::vector<size_t> toremove;
28  for (size_t i = 0; i < N; ++i)
29  dx[i] = fastabs(toadd[i]);
30  // First clean up new vectors
31  size_t first = 0;
32  for (size_t i = 1; i < N; ++i)
33  {
34  // Are they the same index?
35  if (duplicate[i] == duplicate[first])
36  {
37  for (size_t j = first; j < i; ++j)
38  {
39  //Are they the same direction
40  double dxmax = std::max(std::max(dx[first], dx[j]), 1e-30);
41  if (fastabs(toadd[j] - toadd[i]) < dxmax*0.1)
42  {
43  toremove.push_back(i);
44  break;
45  }
46  }
47  }
48  else
49  first = i;
50  }
51  RemoveVector(duplicate, toremove);
52  RemoveVector(toadd, toremove);
53  }
54 
55  void PeriodicGetRidDuplicates(std::vector<int> &duplicate, std::vector<Vector2D> &toadd,
56  std::vector<int> const& old_duplicate, std::vector<Vector2D> const& old_toadd,double Rproc)
57  {
58  if (duplicate.empty())
59  return;
60  assert(duplicate.size() == toadd.size());
61  std::vector<size_t> indeces = sort_index(duplicate);
62  toadd = VectorValues(toadd, indeces);
63  duplicate = VectorValues(duplicate, indeces);
64  size_t N = duplicate.size();
65  size_t Nold = old_duplicate.size();
66  std::vector<double> dx(N, 0), dx_old(Nold, 0);
67  std::vector<size_t> toremove;
68  for (size_t i = 0; i < N; ++i)
69  dx[i] = fastabs(toadd[i]);
70  for (size_t i = 0; i < Nold; ++i)
71  dx_old[i] = fastabs(old_toadd[i]);
72  // First clean up new vectors
73  PeriodicGetRidDuplicatesSingle(duplicate, toadd);
74  // Now clean up with respect to old vectors
75  toremove.clear();
76  N = duplicate.size();
77  for (size_t i = 0; i < N; ++i)
78  {
79  // Are they the same index?
80  size_t index = static_cast<size_t>(std::lower_bound(old_duplicate.begin(), old_duplicate.end(), duplicate[i])
81  - old_duplicate.begin());
82  if (index < Nold)
83  {
84  for (size_t j = index; j < Nold; ++j)
85  {
86  if (old_duplicate[j] != duplicate[i])
87  break;
88  //Are they the same direction
89  double dxmax = std::max(std::max(dx[i], dx_old[j]), Rproc*0.1);
90  if (fastabs(old_toadd[j] - toadd[i]) < dxmax*0.1)
91  {
92  toremove.push_back(i);
93  break;
94  }
95  }
96  }
97  }
98  RemoveVector(duplicate, toremove);
99  RemoveVector(toadd, toremove);
100  }
101 
102  template<class T> bool is_in
103  (const T& t,
104  const vector<T>& v)
105  {
106  BOOST_FOREACH(const T&m, v) {
107  if (t == m)
108  return true;
109  }
110  return false;
111  }
112 #endif
113 
114  pair<int, int> find_diff(const facet& f1, const facet& f2)
115  {
116  if (f1.vertices.first != f2.vertices.first &&
117  f1.vertices.first != f2.vertices.second &&
118  f1.vertices.first != f2.vertices.third)
119  return pair<int, int>(f1.vertices.first, 0);
120  else if (f1.vertices.second != f2.vertices.first &&
121  f1.vertices.second != f2.vertices.second &&
122  f1.vertices.second != f2.vertices.third)
123  return pair<int, int>(f1.vertices.second, 1);
124  else if (f1.vertices.third != f2.vertices.first &&
125  f1.vertices.third != f2.vertices.second &&
126  f1.vertices.third != f2.vertices.third)
127  return pair<int, int>(f1.vertices.third, 2);
128  else
129  throw UniversalError("Delaunay, Couldn't find difference bewteen two facets");
130  }
131 }
132 
133 Delaunay::DataOnlyForBuild::DataOnlyForBuild() :copied(vector<vector<char> >())
134 {}
135 
136 Delaunay::DataOnlyForBuild::DataOnlyForBuild(DataOnlyForBuild const& other) : copied(other.copied) {}
137 
138 Delaunay::DataOnlyForBuild& Delaunay::DataOnlyForBuild::operator=
139 (DataOnlyForBuild const& other)
140 {
141  if (this != &other)
142  {
143  copied = other.copied;
144  }
145  return *this;
146 }
147 
148 Delaunay::Delaunay(void) :
149  lastFacet(0), CalcRadius(false),
150  radius(vector<double>()),
151  PointWasAdded(false),
152  last_facet_added(0),
153  f(vector<facet>()),
154  cor(vector<Vector2D>()),
155  length(0),
156  olength(0), location_pointer(0), last_loc(0),
157  OrgIndex(vector<int>()),
158  logger(0)
159 {}
160 
161 Delaunay::Delaunay(Delaunay const& other) :
162  lastFacet(other.lastFacet),
163  CalcRadius(other.CalcRadius),
164  radius(other.radius),
165  PointWasAdded(other.PointWasAdded),
166  last_facet_added(other.last_facet_added),
167  f(other.f),
168  cor(other.cor),
169  length(other.length),
170  olength(other.olength),
171  location_pointer(other.location_pointer),
172  last_loc(other.last_loc),
173  OrgIndex(other.OrgIndex),
174  logger(other.logger)
175 {}
176 
177 Delaunay::~Delaunay(void)
178 {
179  cor.clear();
180  f.clear();
181 }
182 
183 namespace
184 {
185  // Checks if a point is inside a triangle
186  bool InTriangle(const TripleConstRef<Vector2D>& tri,
187  const Vector2D& point)
188  {
189  return (orient2d(TripleConstRef<Vector2D>(tri.first,
190  tri.second,
191  point)) > 0) &&
192  (orient2d(TripleConstRef<Vector2D>(tri.second,
193  tri.third,
194  point)) > 0) &&
195  (orient2d(TripleConstRef<Vector2D>(tri.third,
196  tri.first,
197  point)) > 0);
198  }
199 
200  vector<double> CellSize(std::vector<std::pair<Vector2D, Vector2D> > const& points)
201  {
202  int n = static_cast<int>(points.size());
203  double minx = points[0].first.x;
204  double miny = points[0].first.y;
205  double maxx = minx;
206  double maxy = miny;
207  for (int i = 1; i < n; ++i)
208  {
209  minx = min(points[static_cast<size_t>(i)].first.x, minx);
210  miny = min(points[static_cast<size_t>(i)].first.y, miny);
211  maxx = max(points[static_cast<size_t>(i)].first.x, maxx);
212  maxy = max(points[static_cast<size_t>(i)].first.y, maxy);
213  }
214  vector<double> res(4);
215  res[0] = minx;
216  res[1] = maxx;
217  res[2] = miny;
218  res[3] = maxy;
219  return res;
220  }
221 
222  size_t find_index(facet const& fc, int i)
223  {
224  for (size_t j = 0; j < 3; ++j)
225  {
226  if (fc.neighbors[j] == i)
227  return j;
228  }
229  throw UniversalError("Error in find_index: Index not found");
230  }
231 }
232 
233 void Delaunay::add_point(size_t index, stack<std::pair<size_t, size_t> > &flip_stack)
234 {
235  // Check if point is inside big triangle
236 
237  const size_t triangle = Walk(index);
238  const Triplet<int> outer(f[triangle].vertices);
239  const Triplet<int> temp_friends(f[triangle].neighbors);
240  f[triangle].vertices.set(outer.third, outer.first, static_cast<int>(index));
241  f[triangle].neighbors.set(temp_friends.third,
242  location_pointer + 1,
243  location_pointer + 2);
244  f.push_back(facet(TripleConstRef<int>(outer.first,
245  outer.second,
246  static_cast<int>(index)),
247  TripleConstRef<int>(temp_friends.first,
248  location_pointer + 2,
249  static_cast<int>(triangle))));
250  f.push_back(facet(TripleConstRef<int>(outer.second,
251  outer.third,
252  static_cast<int>(index)),
253  TripleConstRef<int>(temp_friends.second,
254  static_cast<int>(triangle),
255  location_pointer + 1)));
256  // _update the friends list of the friends
257  if (temp_friends.second != last_loc)
258  {
259  const size_t i = find_index(f[static_cast<size_t>(temp_friends.second)], static_cast<int>(triangle));
260  f[static_cast<size_t>(temp_friends.second)].neighbors[i] = location_pointer + 2;
261  }
262  if (temp_friends.first != last_loc)
263  {
264  const size_t i = find_index(f[static_cast<size_t>(temp_friends.first)], static_cast<int>(triangle));
265  f[static_cast<size_t>(temp_friends.first)].neighbors[i] = location_pointer + 1;
266  }
267  // Calculate radius if needed
268  if (CalcRadius)
269  {
270  radius[static_cast<size_t>(triangle)] = CalculateRadius(static_cast<int>(triangle));
271  int n = int(f.size());
272  int m = int(radius.size());
273  if (n > m - 1)
274  {
275  radius.push_back(CalculateRadius(location_pointer + 1));
276  radius.push_back(CalculateRadius(location_pointer + 2));
277  }
278  else
279  if (n > m)
280  {
281  radius[static_cast<size_t>(location_pointer) + 1] = CalculateRadius(location_pointer + 1);
282  radius.push_back(CalculateRadius(location_pointer + 2));
283  }
284  else
285  {
286  radius[static_cast<size_t>(location_pointer) + 1] = CalculateRadius(location_pointer + 1);
287  radius[static_cast<size_t>(location_pointer) + 2] = CalculateRadius(location_pointer + 2);
288  }
289  }
290 
291  // check if flipping is needed
292  flip(triangle, static_cast<size_t>(temp_friends.third), flip_stack);
293  flip(static_cast<size_t>(location_pointer) + 1, static_cast<size_t>(temp_friends.first), flip_stack);
294  flip(static_cast<size_t>(location_pointer) + 2, static_cast<size_t>(temp_friends.second), flip_stack);
295 
296  // _update number of facets
297  location_pointer += 2;
298 }
299 
300 void Delaunay::flip(size_t i, size_t j, stack<std::pair<size_t, size_t> > &flip_stack)
301 {
302  if (j == static_cast<size_t>(last_loc))
303  return;
304  flip_stack.push(std::pair<size_t, size_t>(i, j));
305  while (!flip_stack.empty())
306  {
307  const pair<size_t, size_t> indexes = flip_stack.top();
308  // Returns the index to the point to check in coordinates and the index of the point in the facet
309  const pair<int, int> check = find_diff(f[indexes.second],
310  f[indexes.first]);
311  const pair<int, int> other = find_diff(f[indexes.first],
312  f[indexes.second]);
313 
314  facet& prefetch_1 = f[indexes.first];
315  if (incircle(cor[static_cast<size_t>(prefetch_1.vertices.first)],
316  cor[static_cast<size_t>(prefetch_1.vertices.second)],
317  cor[static_cast<size_t>(prefetch_1.vertices.third)],
318  cor[static_cast<size_t>(check.first)]) > 0)
319  {
320  //The point is in a circle change the facets and their friends
321  const int v1 = prefetch_1.vertices[static_cast<size_t>(other.second + 1) % 3];
322  const int f1 = prefetch_1.neighbors[static_cast<size_t>(other.second)];
323  const int f12 = prefetch_1.neighbors[static_cast<size_t>(other.second + 2) % 3];
324  facet& prefetch_2 = f[indexes.second];
325  const int v2 = prefetch_2.vertices[static_cast<size_t>(check.second + 1) % 3];
326  const int f2 = prefetch_2.neighbors[static_cast<size_t>(check.second + 2) % 3];
327  const int f22 = prefetch_2.neighbors[static_cast<size_t>(check.second)];
328  prefetch_1.vertices.set(other.first, v1, check.first);
329  prefetch_2.vertices.set(check.first, v2, other.first);
330  prefetch_1.neighbors.set(f1, f2, static_cast<int>(indexes.second));
331  prefetch_2.neighbors.set(f22, f12, static_cast<int>(indexes.first));
332  // change the friends of the friends if needed
333  if (f2 != last_loc)
334  {
335  f[static_cast<size_t>(f2)].neighbors[static_cast<size_t>(find_index(f[static_cast<size_t>(f2)], static_cast<int>(indexes.second)))] = static_cast<int>(indexes.first);
336  }
337  if (f12 != last_loc)
338  {
339  f[static_cast<size_t>(f12)].neighbors[static_cast<size_t>(find_index(f[static_cast<size_t>(f12)], static_cast<int>(indexes.first)))] = static_cast<int>(indexes.second);
340  }
341  // Calculate the new radius if needed
342  if (CalcRadius)
343  {
344  radius[indexes.first] = CalculateRadius(static_cast<int>(indexes.first));
345  radius[indexes.second] = CalculateRadius(static_cast<int>(indexes.second));
346  }
347  // clear the checked facets
348  flip_stack.pop();
349  // push into the stack the new facets to check
350  if (prefetch_2.neighbors.first != last_loc)
351  flip_stack.push(std::pair<size_t, size_t>(indexes.second,
352  static_cast<size_t>(prefetch_2.neighbors.first)));
353  if (prefetch_1.neighbors.second != last_loc)
354  flip_stack.push(std::pair<size_t, size_t>(indexes.first,
355  static_cast<size_t>(prefetch_1.neighbors.second)));
356  }
357  else
358  {
359  // clear the checked facets
360  flip_stack.pop();
361  }
362  }
363 }
364 
365 namespace{
366  /*
367  int get_rank(void)
368  {
369  int rank = -1;
370  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
371  return rank;
372  }
373  */
374 }
375 
376 void Delaunay::build_delaunay(vector<Vector2D>const& vp, std::vector<std::pair<Vector2D, Vector2D> > const& cpoints)
377 {
378  DataOnlyForBuild data;
379  lastFacet = 0;
380  CalcRadius = false;
381  length = int(vp.size() + 3);
382  int len = length - 3;
383  olength = static_cast<size_t>(len);
384  f.clear();
385  cor.clear();
386  f.reserve(static_cast<size_t>(2 * length + 1 + static_cast<int>(17 * sqrt(1.*length))));
387  cor.reserve(static_cast<size_t>(length + 9 * static_cast<int>(sqrt(1.*length))));
388  last_loc = INT_MAX;
389  for (int i = 0; i < len; i++)
390  {
391  cor.push_back(vp[static_cast<size_t>(i)]);
392  }
393 
394  // add the 3 extreme points
395  Vector2D p_temp;
396  vector<double> cellsize = CellSize(cpoints);
397  double width = cellsize[1] - cellsize[0];
398  double height = cellsize[3] - cellsize[2];
399  width = max(width, height);
400  height = max(width, height);
401  p_temp.x = cellsize[0] - 100 * width;
402  p_temp.y = cellsize[2] - 100 * height;
403  cor.push_back(p_temp);
404  p_temp.x = cellsize[1] + 100 * width;
405  p_temp.y = cellsize[2] - 100 * height;
406  cor.push_back(p_temp);
407  p_temp.x = (cellsize[0] + cellsize[1]) / 2.0;
408  p_temp.y = cellsize[3] + 100 * height;
409  cor.push_back(p_temp);
410  // Create the big triangle, and assign friends
411  facet f_temp;
412  f.push_back(f_temp);
413  f[0].vertices[0] = len;
414  f[0].vertices[1] = len + 1;
415  f[0].vertices[2] = len + 2;
416  for (size_t i = 0; i < 3; i++)
417  f[0].neighbors[i] = last_loc;
418  location_pointer = 0;
419  // add the points
420  size_t nloop = static_cast<size_t>(length) - 3;
421  stack<std::pair<size_t, size_t> > flip_stack;
422  for (size_t i = 0; i < nloop; i++)
423  add_point(i, flip_stack);
424  // Calculate radius
425  radius.resize(f.size());
426  int n = int(f.size());
427  for (int i = 0; i < n; ++i)
428  radius[static_cast<size_t>(i)] = CalculateRadius(i);
429  //CalcRadius = true;
430 
431 }
432 
433 bool Delaunay::CheckCorrect(void)
434 {
435  std::size_t Ntri = f.size();
436  for (std::size_t i = 0; i < Ntri; ++i)
437  {
438  facet const& T = f[i];
439  for (std::size_t j = 0; j < 3; ++j)
440  {
441  // Check same neighbors
442  if (T.neighbors[j] != last_loc)
443  {
444  bool found = false;
445  for (size_t k = 0; k < 3; ++k)
446  if (f[static_cast<size_t>(T.neighbors[j])].neighbors[k] == static_cast<int>(i))
447  found = true;
448  if (!found)
449  {
450  UniversalError eo("Failed checkcorrect in delaunay");
451  eo.AddEntry("Facet checked", static_cast<double>(i));
452  eo.AddEntry("Neighbor that is not dual neighbor", static_cast<double>(T.neighbors[j]));
453  throw eo;
454  }
455  }
456  // Check incircle
457  if (T.neighbors[j] != last_loc)
458  {
459  // Find point to shared
460  for (size_t k = 0; k < 3; ++k)
461  {
462  bool found = false;
463  for (size_t l = 0; l < 3; ++l)
464  {
465  if (T.vertices[l] == f[static_cast<size_t>(T.neighbors[j])].vertices[k])
466  {
467  found = true;
468  break;
469  }
470  }
471  if (!found)
472  {
473  if (incircle(cor[static_cast<size_t>(T.vertices[0])], cor[static_cast<size_t>(T.vertices[1])],
474  cor[static_cast<size_t>(T.vertices[2])], cor[static_cast<size_t>(f[static_cast<size_t>(T.neighbors[j])].vertices[k])]) > 0)
475  {
476  UniversalError eo("Failed checkcorrect in delaunay");
477  eo.AddEntry("Facet checked", static_cast<double>(i));
478  eo.AddEntry("Point that is inside triangle", static_cast<double>(f[static_cast<size_t>(T.neighbors[j])].vertices[k]));
479  throw eo;
480  }
481  break;
482  }
483  }
484  }
485  }
486  }
487  return true;
488 }
489 
490 
491 double Delaunay::triangle_area(int index)
492 {
493  const TripleConstRef<Vector2D> p
494  (cor[static_cast<size_t>(f[static_cast<size_t>(index)].vertices.first)],
495  cor[static_cast<size_t>(f[static_cast<size_t>(index)].vertices.second)],
496  cor[static_cast<size_t>(f[static_cast<size_t>(index)].vertices.third)]);
497  const double x1 = p.third.x - p.first.x;
498  const double x2 = p.second.x - p.first.x;
499  const double y1 = p.third.y - p.first.y;
500  const double y2 = p.second.y - p.first.y;
501  return -0.5*(x1*y2 - x2*y1);
502 }
503 
504 void Delaunay::update(const vector<Vector2D>& points, std::vector<std::pair<Vector2D, Vector2D> >const& cpoints)
505 {
506  if (logger)
507  logger->output(cor, f);
508  build_delaunay(points, cpoints);
509 }
510 
511 namespace {
512 
513  int Triplet<int>::* walk_condition(const vector<Vector2D>& cor,
514  const Triplet<int>& vertices,
515  size_t point)
516  {
517  if (orient2d(TripleConstRef<Vector2D>
518  (cor[static_cast<size_t>(vertices.first)],
519  cor[static_cast<size_t>(vertices.second)],
520  cor[point])) < 0)
521  return &Triplet<int>::first;
522  else if (orient2d(TripleConstRef<Vector2D>
523  (cor[static_cast<size_t>(vertices.second)],
524  cor[static_cast<size_t>(vertices.third)],
525  cor[point])) < 0)
526  return &Triplet<int>::second;
527  else if (orient2d(TripleConstRef<Vector2D>
528  (cor[static_cast<size_t>(vertices.third)],
529  cor[static_cast<size_t>(vertices.first)],
530  cor[point])) < 0)
531  return &Triplet<int>::third;
532  else
533  return 0;
534  }
535 
536  class WalkBookkeeper
537  {
538  public:
539 
540  private:
541 
542  };
543 
544  size_t find_new_facet(const vector<Vector2D>& cor,
545  const vector<facet>& f,
546  size_t point,
547  size_t last_facet)
548  {
549  size_t res = last_facet;
550  int Triplet<int>::* next = walk_condition(cor,
551  f[res].vertices,
552  point);
553  while (next) {
554  res = static_cast<size_t>(f[res].neighbors.*next);
555  next = walk_condition(cor,
556  f[res].vertices,
557  point);
558  }
559  return res;
560  }
561 }
562 
563 size_t Delaunay::Walk(size_t point)
564 {
565  lastFacet = static_cast<int>(find_new_facet(cor, f, point, static_cast<size_t>(lastFacet)));
566  return static_cast<size_t>(lastFacet);
567 }
568 
569 vector<int> Delaunay::FindContainingTetras(int StartTetra, int point)
570 {
571  vector<int> res;
572  FindContainingTetras(StartTetra, point, res);
573  return res;
574 }
575 
576 double Delaunay::FindMaxRadius(int point)
577 {
578  const vector<int> vec = FindContainingTetras(static_cast<int>(Walk(static_cast<size_t>(point))), point);
579  double r = 0;
580  for (size_t i = 0; i < vec.size(); ++i)
581  r = max(r, radius[static_cast<size_t>(vec[static_cast<size_t>(i)])]);
582  return 2 * r;
583 }
584 
585 void Delaunay::FindContainingTetras(int StartFacet, int point, vector<int> &result)
586 {
587  result.clear();
588  int PointLocation = FindPointInFacet(StartFacet, point);
589  int NextFacet = f[static_cast<size_t>(StartFacet)].neighbors[static_cast<size_t>(PointLocation)];
590  result.reserve(12);
591  result.push_back(NextFacet);
592  while (NextFacet != StartFacet)
593  {
594  PointLocation = FindPointInFacet(NextFacet, point);
595  NextFacet = f[static_cast<size_t>(NextFacet)].neighbors[static_cast<size_t>(PointLocation)];
596  result.push_back(NextFacet);
597  }
598 }
599 
600 int Delaunay::FindPointInFacet(int facet, int point)
601 {
602  for (int i = 0; i < 3; ++i)
603  if (f[static_cast<size_t>(facet)].vertices[static_cast<size_t>(i)] == point)
604  return i;
605  UniversalError eo("Error in Delaunay, FindPointInFacet");
606  eo.AddEntry("Facet number", facet);
607  eo.AddEntry("Point number", point);
608  throw eo;
609 }
610 
611 bool Delaunay::IsOuterFacet(int facet)const
612 {
613  //int PointNum=length-1;
614  for (int i = 0; i < 3; ++i)
615  for (size_t j = 0; j < 3; ++j)
616  if (f[static_cast<size_t>(facet)].vertices[static_cast<size_t>(i)] == static_cast<int>(olength + j))
617  return true;
618  return false;
619 }
620 
621 double Delaunay::CalculateRadius(int facet)
622 {
623  /*const double small = 1e-8;
624  const double a = cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[0])].distance(cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[1])]);
625  const double b = cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[0])].distance(cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[2])]);
626  const double c = cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[2])].distance(cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[1])]);
627  const double temp1 = b + c - a;
628  const double temp2 = c + a - b;
629  const double temp3 = b - c + a;
630  double maxv = std::max(std::max(std::abs(temp1), std::abs(temp2)), std::abs(temp3));
631  if (std::abs(temp1) < (maxv*small) || std::abs(temp2) < (maxv*small) || std::abs(temp3) < (maxv*small))
632  return CalcRadiusHiRes(facet);
633  return a*b*c / std::sqrt((a + b + c)*temp1*temp2*temp3);*/
634  double R = 0;
635  GetCircleCenter(facet, R);
636  return R;
637 }
638 
639 double Delaunay::CalcRadiusHiRes(int facet)
640 {
641  std::array<boost::multiprecision::cpp_dec_float_50, 2> V0,V1,V2;
642  V0[0] = cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[0])].x;
643  V0[1] = cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[0])].y;
644  V1[0] = cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[1])].x;
645  V1[1] = cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[1])].y;
646  V2[0] = cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[2])].x;
647  V2[1] = cor[static_cast<size_t>(f[static_cast<size_t>(facet)].vertices[2])].y;
648 
649  boost::multiprecision::cpp_dec_float_50 a = boost::multiprecision::sqrt((V0[0] - V1[0])*(V0[0] - V1[0]) + (V0[1] - V1[1])*(V0[1] - V1[1]));
650  boost::multiprecision::cpp_dec_float_50 b = boost::multiprecision::sqrt((V0[0] - V2[0])*(V0[0] - V2[0]) + (V0[1] - V2[1])*(V0[1] - V2[1]));
651  boost::multiprecision::cpp_dec_float_50 c = boost::multiprecision::sqrt((V2[0] - V1[0])*(V2[0] - V1[0]) + (V2[1] - V1[1])*(V2[1] - V1[1]));
652 
653  boost::multiprecision::cpp_dec_float_50 temp1 = b + c - a;
654  boost::multiprecision::cpp_dec_float_50 temp2 = -b + c + a;
655  boost::multiprecision::cpp_dec_float_50 temp3 = b - c - a;
656  boost::multiprecision::cpp_dec_float_50 temp4 = a + b + c;
657  temp4 = a * b*c / boost::multiprecision::sqrt(temp4*temp1*temp2*temp3);
658  return (temp4).convert_to<double>();
659 }
660 
661 int Delaunay::GetOriginalIndex(int NewPoint) const
662 {
663  return NewPoint;
664 }
665 
666 double Delaunay::GetFacetRadius(int facet) const
667 {
668  return radius[static_cast<size_t>(facet)];
669 }
670 
671 void Delaunay::ChangeOlength(int n)
672 {
673  olength = static_cast<size_t>(n);
674 }
675 
676 void Delaunay::Changelength(int n)
677 {
678  length = n + 3;
679 }
680 
681 vector<Vector2D>& Delaunay::ChangeCor(void)
682 {
683  return cor;
684 }
685 
686 const vector<Vector2D>& Delaunay::getCor(void) const
687 {
688  return cor;
689 }
690 
691 const facet& Delaunay::get_facet(int index) const
692 {
693  return f[static_cast<size_t>(index)];
694 }
695 
696 double Delaunay::get_facet_coordinate(int Facet, int vertice, int dim)
697 {
698  if (dim == 0)
699  return cor[static_cast<size_t>(f[static_cast<size_t>(Facet)].vertices[static_cast<size_t>(vertice)])].x;
700  else
701  return cor[static_cast<size_t>(f[static_cast<size_t>(Facet)].vertices[static_cast<size_t>(vertice)])].y;
702 }
703 
704 Vector2D Delaunay::get_point(size_t index) const
705 {
706  return cor[index];
707 }
708 
709 double Delaunay::get_cor(int index, int dim) const
710 {
711  if (dim == 0)
712  return cor[static_cast<size_t>(index)].x;
713  else if (dim == 1)
714  return cor[static_cast<size_t>(index)].y;
715  else
716  throw UniversalError("Error in Delaunay::get_cor. Invalid index");
717 }
718 
719 int Delaunay::get_num_facet(void) const
720 {
721  return static_cast<int>(f.size());
722 }
723 
724 int Delaunay::get_length(void) const
725 {
726  return length - 3;
727 }
728 
729 int Delaunay::get_last_loc(void) const
730 {
731  return last_loc;
732 }
733 
734 void Delaunay::set_point(int index, Vector2D p)
735 {
736  cor[static_cast<size_t>(index)] = p;
737 }
738 
739 int Delaunay::GetOriginalLength(void) const
740 {
741  return static_cast<int>(olength);
742 }
743 
744 vector<Vector2D>& Delaunay::GetMeshPoints(void)
745 {
746  return cor;
747 }
748 
749 int Delaunay::GetTotalLength(void)
750 {
751  return static_cast<int>(cor.size());
752 }
753 
754 void Delaunay::AddBoundaryPoints(vector<Vector2D> const& points)
755 {
756  int n = static_cast<int>(points.size());
757  int N = static_cast<int>(cor.size());
758  stack<std::pair<size_t, size_t> > flip_stack;
759 // /*vector<int> order=*/HilbertOrder(points,n);
760  cor.insert(cor.end(), points.begin(), points.end());
761  TripleConstRef<Vector2D> OuterTri(cor[olength], cor[olength + 1], cor[olength + 2]);
762  for (int i = 0; i < n; ++i)
763  if(InTriangle(OuterTri,cor[static_cast<size_t>(N + i)]))
764  add_point(static_cast<size_t>(N + i), flip_stack);
765 }
766 
767 void Delaunay::AddAditionalPoint(Vector2D const& vec)
768 {
769  cor.push_back(vec);
770 }
771 
772 int Delaunay::GetCorSize(void)const
773 {
774  return static_cast<int>(cor.size());
775 }
776 
777 bool Delaunay::IsTripleOut(int index) const
778 {
779  int counter = 0;
780  for (size_t i = 0; i < 3; ++i)
781  if (IsOuterFacet(f[static_cast<size_t>(index)].neighbors[static_cast<size_t>(i)]))
782  ++counter;
783  if (counter > 1)
784  return true;
785  else
786  return false;
787 }
788 
789 int Delaunay::FindTripleLoc(facet const& fct)const
790 {
791  for (size_t i = 0; i < 3; ++i)
792  if (!IsOuterFacet(fct.neighbors[static_cast<size_t>(i)]))
793  return static_cast<int>((i + 1) % 3);
794  throw UniversalError("Trouble in constructing boundary triangles. No inner neighbor");
795 }
796 
797 namespace
798 {
799  bool IsOuterQuick(facet const& f, int olength)
800  {
801  for (size_t i = 0; i < 3; ++i)
802  if (f.vertices[static_cast<size_t>(i)] >= olength)
803  return true;
804  return false;
805  }
806 
807  bool IsEdgeFacet(vector<facet> const& facets, facet const& f, int olength)
808  {
809  int counter = 0;
810  for (size_t i = 0; i < 3; ++i)
811  {
812  if (f.vertices[static_cast<size_t>(i)] >= olength)
813  return false;
814  if (IsOuterQuick(facets[static_cast<size_t>(f.neighbors[static_cast<size_t>(i)])], olength))
815  ++counter;
816  }
817  if (counter > 0)
818  return true;
819  else
820  return false;
821  }
822 
823  bool CircleSegmentIntersect(Edge const& edge, Vector2D const& center, double R)
824  {
825  Vector2D AC = center - edge.vertices.first;
826  Vector2D AB = edge.vertices.second - edge.vertices.first;
827  double d = ScalarProd(AC, AB);
828  if (d < 0)
829  {
830  if (abs(AC) > R)
831  return false;
832  else
833  return true;
834  }
835  double LAB = abs(AB);
836  if (d > LAB*LAB)
837  {
838  if (abs(center - edge.vertices.second) > R)
839  return false;
840  else
841  return true;
842  }
843  Vector2D closest = edge.vertices.first + AB*d / (LAB*LAB);
844  if (abs(center - closest) > R)
845  return false;
846  else
847  return true;
848  }
849 }
850 
851 vector<int> Delaunay::GetOuterFacets(int start_facet, int real_point, int olength2)
852 {
853  int cur_facet = start_facet;
854  vector<int> f_temp, containing_facets;
855  f_temp.reserve(static_cast<size_t>(10 * sqrt(1.0*olength2)));
856  int point_index = FindPointInFacet(cur_facet, real_point);
857  if (IsOuterQuick(f[static_cast<size_t>(f[static_cast<size_t>(cur_facet)].neighbors[static_cast<size_t>(point_index)])], olength2))
858  {
859  point_index = (point_index + 1) % 3;
860  real_point = f[static_cast<size_t>(cur_facet)].vertices[static_cast<size_t>(point_index)];
861  }
862  if (IsOuterQuick(f[static_cast<size_t>(f[static_cast<size_t>(cur_facet)].neighbors[static_cast<size_t>(point_index)])], olength2))
863  {
864  point_index = (point_index + 1) % 3;
865  real_point = f[static_cast<size_t>(cur_facet)].vertices[static_cast<size_t>(point_index)];
866  }
867  do
868  {
869  FindContainingTetras(cur_facet, real_point, containing_facets);
870  int old_current = cur_facet;
871  for (size_t i = 0; i < containing_facets.size(); ++i)
872  {
873  if (IsEdgeFacet(f, f[static_cast<size_t>(containing_facets[static_cast<size_t>(i)])], olength2) &&
874  containing_facets[static_cast<size_t>(i)] != old_current)
875  cur_facet = containing_facets[static_cast<size_t>(i)];
876  if (!IsOuterQuick(f[static_cast<size_t>(containing_facets[static_cast<size_t>(i)])], olength2))
877  f_temp.push_back(containing_facets[static_cast<size_t>(i)]);
878  }
879  point_index = (1 + FindPointInFacet(cur_facet, real_point)) % 3;
880  if (IsTripleOut(cur_facet))
881  point_index = (point_index + 1) % 3;
882  real_point = f[static_cast<size_t>(cur_facet)].vertices[static_cast<size_t>(point_index)];
883  } while (start_facet != cur_facet);
884  sort(f_temp.begin(), f_temp.end());
885  f_temp = unique(f_temp);
886  return f_temp;
887 }
888 
889 namespace {
890 
891  vector<int> range(int len)
892  {
893  vector<int> res(static_cast<size_t>(len));
894  for(int i=0;i<len;++i)
895  res.at(static_cast<size_t>(i)) = i;
896  return res;
897  }
898 
899  vector<vector<int> > replicate_all_points(const vector<Edge>& edges,
900  size_t olength)
901  {
902  vector<vector<int> > res(edges.size());
903  const vector<int> my_range = range(static_cast<int>(olength));
904  for(size_t i=0;i<edges.size(); ++i)
905  res.at(i) = my_range;
906  return res;
907  }
908 }
909 
910  namespace {
911 
912  void clean_toduplicate(vector<vector<int> >& toduplicate)
913  {
914  for(size_t i=0;i<toduplicate.size();++i){
915  sort(toduplicate.at(i).begin(),
916  toduplicate.at(i).end());
917  toduplicate.at(i) = unique(toduplicate.at(i));
918  }
919  }
920  }
921 
922 vector<vector<int> > Delaunay::FindOuterPoints(vector<Edge> const& edges)
923 {
924  if(olength<100)
925  return replicate_all_points(edges, olength);
926 
927  // Walk to an outer point
928  vector<vector<int> > toduplicate(edges.size());
929  {
930  int cur_facet = static_cast<int>(Walk(olength));
931  vector<bool> checked(f.size(), false);
932  AddOuterFacets(cur_facet, toduplicate, edges, checked);
933  }
934  clean_toduplicate(toduplicate);
935  return toduplicate;
936 }
937 
938 void Delaunay::AddPeriodicMPI(std::vector<int> &toduplicate, std::vector<Vector2D> &periodic_add)
939 {
940  assert(periodic_add.size() == toduplicate.size());
941  std::vector<int> toremove;
942  std::vector<Vector2D> toadd;
943  for (size_t j = 0; j < toduplicate.size(); ++j)
944  {
945  Vector2D temp = cor[static_cast<size_t>(toduplicate[j])] + periodic_add[j];
946  if (InTriangle(TripleConstRef<Vector2D>
947  (cor[static_cast<size_t>(olength)],
948  cor[static_cast<size_t>(olength + 1)],
949  cor[static_cast<size_t>(olength + 2)]),
950  temp))
951  toadd.push_back(temp);
952  else
953  toremove.push_back(static_cast<int>(j));
954  }
955  RemoveVector(toduplicate, toremove);
956  try
957  {
958  AddBoundaryPoints(toadd);
959  }
960  catch (UniversalError &eo)
961  {
962  eo.AddEntry("Error in AddRigid", 0);
963  throw;
964  }
965 }
966 
967 void Delaunay::AddRigid(vector<Edge> const& edges,
968  vector<vector<int> > &toduplicate)
969 {
970  vector<int> toremove;
971  for (size_t i = 0; i < edges.size(); ++i)
972  {
973  toremove.clear();
974  if (toduplicate[i].empty())
975  continue;
976  vector<Vector2D> toadd;
977  toadd.reserve(toduplicate[i].size());
978  Vector2D par(Parallel(edges[i]));
979  par = par / abs(par);
980  for (size_t j = 0; j < toduplicate[i].size(); ++j)
981  {
982  Vector2D temp = cor[static_cast<size_t>(toduplicate[i][j])] - edges[i].vertices.first;
983  temp = 2 * par*ScalarProd(par, temp) - temp + edges[i].vertices.first;
984  if (InTriangle(TripleConstRef<Vector2D>
985  (cor[static_cast<size_t>(olength)],
986  cor[static_cast<size_t>(olength + 1)],
987  cor[static_cast<size_t>(olength + 2)]),
988  temp))
989  toadd.push_back(temp);
990  else
991  toremove.push_back(static_cast<int>(j));
992  }
993  RemoveVector(toduplicate[i], toremove);
994  try
995  {
996  AddBoundaryPoints(toadd);
997  }
998  catch (UniversalError &eo)
999  {
1000  eo.AddEntry("Error in AddRigid", 0);
1001  throw;
1002  }
1003  }
1004 }
1005 
1006 namespace
1007 {
1008  vector<Edge> GetCornerEdges(OuterBoundary const* obc)
1009  {
1010  const double dx = obc->GetGridBoundary(Right) - obc->GetGridBoundary(Left);
1011  const double dy = obc->GetGridBoundary(Up) - obc->GetGridBoundary(Down);
1012  vector<Edge> res;
1013  const Vector2D RU(obc->GetGridBoundary(Right), obc->GetGridBoundary(Up));
1014  const Vector2D LU(obc->GetGridBoundary(Left), obc->GetGridBoundary(Up));
1015  const Vector2D LD(obc->GetGridBoundary(Left), obc->GetGridBoundary(Down));
1016  const Vector2D RD(obc->GetGridBoundary(Right), obc->GetGridBoundary(Down));
1017  res.push_back(Edge(RU, Vector2D(dx, 0) + RU, 0, 0));
1018  res.push_back(Edge(RU, Vector2D(0, dy) + RU, 0, 0));
1019  res.push_back(Edge(LU, Vector2D(0, dy) + LU, 0, 0));
1020  res.push_back(Edge(LU, Vector2D(-dx, 0) + LU, 0, 0));
1021  res.push_back(Edge(LD, Vector2D(-dx, 0) + LD, 0, 0));
1022  res.push_back(Edge(LD, Vector2D(0, -dy) + LD, 0, 0));
1023  res.push_back(Edge(RD, Vector2D(0, -dy) + RD, 0, 0));
1024  res.push_back(Edge(RD, Vector2D(dx, 0) + RD, 0, 0));
1025  return res;
1026  }
1027 }
1028 
1029 vector<vector<int> > Delaunay::AddPeriodic(OuterBoundary const* obc, vector<Edge> const& edges,
1030  vector<vector<int> > &toduplicate)
1031 {
1032  const double dx = obc->GetGridBoundary(Right) - obc->GetGridBoundary(Left);
1033  const double dy = obc->GetGridBoundary(Up) - obc->GetGridBoundary(Down);
1034  for (size_t i = 0; i < edges.size(); ++i)
1035  {
1036  if (toduplicate[static_cast<size_t>(i)].empty())
1037  continue;
1038  Vector2D change;
1039  switch (i)
1040  {
1041  case(0):
1042  change.x = -dx;
1043  break;
1044  case(1):
1045  change.y = -dy;
1046  break;
1047  case(2):
1048  change.x = dx;
1049  break;
1050  case(3):
1051  change.y = dy;
1052  break;
1053  }
1054  vector<Vector2D> toadd;
1055  toadd.reserve(toduplicate[static_cast<size_t>(i)].size());
1056  //vector<int> pointstemp(toduplicate[static_cast<size_t>(i)].size());
1057  for (size_t j = 0; j < toduplicate[static_cast<size_t>(i)].size(); ++j)
1058  {
1059  toadd.push_back(cor[static_cast<size_t>(toduplicate[static_cast<size_t>(i)][static_cast<size_t>(j)])] + change);
1060  // pointstemp[j]=j;
1061  }
1062  vector<int> order = HilbertOrder(toadd, static_cast<int>(toadd.size()));
1063  ReArrangeVector(toadd, order);
1064  AddBoundaryPoints(toadd);
1065  ReArrangeVector(toduplicate[static_cast<size_t>(i)], order);
1066  //toduplicate[i]=pointstemp;
1067  }
1068  // Done with sides do corners now
1069  vector<Edge> corneredges = GetCornerEdges(obc);
1070  vector<vector<int> > corners(toduplicate.size());
1071  for (size_t i = 0; i < toduplicate.size(); ++i)
1072  {
1073  for (size_t j = 0; j < toduplicate[static_cast<size_t>(i)].size(); ++j)
1074  {
1075  const int facet_loc = static_cast<int>(Walk(static_cast<size_t>(toduplicate[static_cast<size_t>(i)][static_cast<size_t>(j)])));
1076  const Vector2D center = cor[static_cast<size_t>(toduplicate[static_cast<size_t>(i)][static_cast<size_t>(j)])];
1077  const double R = 2 * GetMaxRadius(toduplicate[static_cast<size_t>(i)][static_cast<size_t>(j)], facet_loc);
1078  if (CircleSegmentIntersect(corneredges[2 * static_cast<size_t>(i)], center, R))
1079  corners[static_cast<size_t>(i)].push_back(toduplicate[static_cast<size_t>(i)][static_cast<size_t>(j)]);
1080  if (CircleSegmentIntersect(corneredges[(2 * static_cast<size_t>(i) + 7) % 8], center, R))
1081  corners[(static_cast<size_t>(i) + 3) % 4].push_back(toduplicate[static_cast<size_t>(i)][static_cast<size_t>(j)]);
1082  }
1083  }
1084  for (size_t i = 0; i < corners.size(); ++i)
1085  {
1086  if (corners[static_cast<size_t>(i)].empty())
1087  continue;
1088  sort(corners[static_cast<size_t>(i)].begin(), corners[static_cast<size_t>(i)].end());
1089  corners[static_cast<size_t>(i)] = unique(corners[static_cast<size_t>(i)]);
1090  Vector2D change;
1091  switch (i)
1092  {
1093  case(0):
1094  change.x = -dx;
1095  change.y = -dy;
1096  break;
1097  case(1):
1098  change.y = -dy;
1099  change.x = dx;
1100  break;
1101  case(2):
1102  change.x = dx;
1103  change.y = dy;
1104  break;
1105  case(3):
1106  change.y = dy;
1107  change.x = -dx;
1108  break;
1109  }
1110  vector<Vector2D> toadd;
1111  toadd.reserve(corners[static_cast<size_t>(i)].size());
1112  // vector<int> pointstemp(corners[i].size());
1113  for (size_t j = 0; j < corners[static_cast<size_t>(i)].size(); ++j)
1114  {
1115  toadd.push_back(cor[static_cast<size_t>(corners[static_cast<size_t>(i)][static_cast<size_t>(j)])] + change);
1116  // pointstemp[j]=j;
1117  }
1118  AddBoundaryPoints(toadd);
1119  }
1120  return corners;
1121 }
1122 
1123 void Delaunay::AddHalfPeriodic(OuterBoundary const* obc, vector<Edge> const& edges,
1124  vector<vector<int> > &toduplicate)
1125 {
1126  const double dx = obc->GetGridBoundary(Right) - obc->GetGridBoundary(Left);
1127  // const double dy=obc->GetGridBoundary(Up)-obc->GetGridBoundary(Down);
1128  for (size_t i = 0; i < edges.size(); ++i)
1129  {
1130  if (toduplicate[static_cast<size_t>(i)].empty())
1131  continue;
1132  Vector2D change;
1133  switch (i)
1134  {
1135  case(0):
1136  change.x = -dx;
1137  break;
1138  case(1):
1139  break;
1140  case(2):
1141  change.x = dx;
1142  break;
1143  case(3):
1144  break;
1145  }
1146  vector<Vector2D> toadd;
1147  toadd.reserve(toduplicate[static_cast<size_t>(i)].size());
1148  //vector<int> pointstemp(toduplicate[static_cast<size_t>(i)].size());
1149  Vector2D par(Parallel(edges[static_cast<size_t>(i)]));
1150  par = par / abs(par);
1151  for (size_t j = 0; j < toduplicate[static_cast<size_t>(i)].size(); ++j)
1152  {
1153  Vector2D temp = cor[static_cast<size_t>(toduplicate[static_cast<size_t>(i)][static_cast<size_t>(j)])];
1154  if (i % 2 == 1)
1155  {
1156  temp -= edges[static_cast<size_t>(i)].vertices.first;
1157  temp = 2 * par*ScalarProd(par, temp) - temp + edges[static_cast<size_t>(i)].vertices.first;
1158  }
1159  toadd.push_back(temp + change);
1160  //pointstemp[j]=j;
1161  }
1162  AddBoundaryPoints(toadd);
1163  //toduplicate[i]=pointstemp;
1164  }
1165 }
1166 
1167 void Delaunay::BuildBoundary(OuterBoundary const* obc, vector<Edge> const& edges)
1168 {
1169  vector<vector<int> > toduplicate = FindOuterPoints(edges);
1170  OrgIndex.clear();
1171  if (obc->GetBoundaryType() == Rectengular)
1172  {
1173  AddRigid(edges, toduplicate);
1174  }
1175  else
1176  {
1177  if (obc->GetBoundaryType() == Periodic)
1178  {
1179  vector<vector<int> > corners = AddPeriodic(obc, edges, toduplicate);
1180  for (size_t i = 0; i < 4; ++i)
1181  toduplicate.push_back(corners[static_cast<size_t>(i)]);
1182 #ifdef RICH_MPI
1183  for (size_t i = 0; i < toduplicate.size(); ++i)
1184  for (size_t j = 0; j < toduplicate[i].size(); ++j)
1185  OrgIndex.push_back(toduplicate[i][j]);
1186 #endif
1187  }
1188  else
1189  {
1190  AddHalfPeriodic(obc, edges, toduplicate);
1191  }
1192  }
1193  for (size_t i = 0; i < toduplicate.size(); ++i)
1194  for (size_t j = 0; j < toduplicate[i].size(); ++j)
1195  OrgIndex.push_back(toduplicate[i][j]);
1196  radius.resize(f.size());
1197  int n = int(f.size());
1198  for (int i = 0; i < n; ++i)
1199  radius[static_cast<size_t>(i)] = CalculateRadius(i);
1200 }
1201 
1202 Vector2D Delaunay::GetCircleCenter(int index, double &R)const
1203 {
1204  Vector2D center;
1205  facet const& F = f[static_cast<size_t>(index)];
1206  double x1 = cor[static_cast<size_t>(F.vertices[0])].x;
1207  double y1 = cor[static_cast<size_t>(F.vertices[0])].y;
1208  double x2 = cor[static_cast<size_t>(F.vertices[1])].x;
1209  double y2 = cor[static_cast<size_t>(F.vertices[1])].y;
1210  double x3 = cor[static_cast<size_t>(F.vertices[2])].x;
1211  double y3 = cor[static_cast<size_t>(F.vertices[2])].y;
1212 
1213  // Do we have a case where two point are very close compared to the third?
1214  double d12 = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2);
1215  double d23 = (x3 - x2)*(x3 - x2) + (y3 - y2)*(y3 - y2);
1216  double d13 = (x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3);
1217  double maxv = std::max(std::max(d12, d23), d13);
1218  double small = 1e-10;
1219  double d_inv = (2 * ((x2-x1)*(y3-y1) - (y2-y1) * (x3-x1)));
1220  if (d12 < (small*maxv) || d23 < (small*maxv) || d13 < (small*maxv) || std::abs(d_inv) <small*maxv)
1221  {
1222  std::array<boost::multiprecision::cpp_dec_float_50, 2> V1,V2, V3;
1223  V1[0] = x1;
1224  V1[1] = y1;
1225  V2[0] = x2;
1226  V2[1] = y2;
1227  V3[0] = x3;
1228  V3[1] = y3;
1229  V2[0] -= V1[0];
1230  V3[0] -= V1[0];
1231  V2[1] -= V1[1];
1232  V3[1] -= V1[1];
1233  boost::multiprecision::cpp_dec_float_50 d = 2 * (V2[0] * V3[1] - V2[1] * V3[0]);
1234  boost::multiprecision::cpp_dec_float_50 r_x = (V3[1] * (V2[1] * V2[1] + V2[0] * V2[0]) - V2[1] * (V3[0] * V3[0] + V3[1] * V3[1])) / d;
1235  boost::multiprecision::cpp_dec_float_50 r_y = (-V3[0] * (V2[1] * V2[1] + V2[0] * V2[0]) + V2[0] * (V3[0] * V3[0] + V3[1] * V3[1])) / d;
1236  R = boost::multiprecision::sqrt(r_x*r_x + r_y * r_y).convert_to<double>();
1237  center.Set(r_x.convert_to<double>() + x1, r_y.convert_to<double>() + y1);
1238  }
1239  else
1240  {
1241  d_inv = 1.0 / d_inv;
1242  x2 -= x1;
1243  x3 -= x1;
1244  y2 -= y1;
1245  y3 -= y1;
1246  double x = (y3*(x2*x2 + y2 * y2) - y2 * (x3*x3 + y3 * y3))*d_inv;
1247  double y = (-x3 * (x2*x2 + y2 * y2) + x2 * (x3*x3 + y3 * y3))*d_inv;
1248  R = std::sqrt(x*x + y * y);
1249  center.Set(x + x1,y + y1);
1250  }
1251  return center;
1252 }
1253 
1254 double Delaunay::GetMaxRadius(int point, int startfacet)
1255 {
1256  double res = 0;
1257  vector<int> neigh = FindContainingTetras(startfacet, point);
1258  for (size_t i = 0; i < neigh.size(); ++i)
1259  res = max(res, radius[static_cast<size_t>(neigh[static_cast<size_t>(i)])]);
1260  return res;
1261 }
1262 
1263 void Delaunay::AddOuterFacets(int tri, vector<vector<int> > &toduplicate,
1264  vector<Edge> const& edges, vector<bool> &checked)
1265 {
1266  stack<int> tocheck;
1267  tocheck.push(tri);
1268  while (!tocheck.empty())
1269  {
1270  int cur_facet = tocheck.top();
1271  tocheck.pop();
1272  for (size_t i = 0; i < 3; ++i)
1273  {
1274  bool added = false;
1275  if (checked[static_cast<size_t>(f[static_cast<size_t>(cur_facet)].vertices[static_cast<size_t>(i)])] || (f[static_cast<size_t>(cur_facet)].vertices[static_cast<size_t>(i)] >= static_cast<int>(olength)))
1276  continue;
1277  vector<int> neigh = FindContainingTetras(cur_facet, f[static_cast<size_t>(cur_facet)].vertices[static_cast<size_t>(i)]);
1278  for (size_t k = 0; k < neigh.size(); ++k)
1279  {
1280  double R = 0;
1281  Vector2D center = GetCircleCenter(neigh[static_cast<size_t>(k)],R);
1282  for (size_t l = 0; l < edges.size(); ++l)
1283  {
1284  if (CircleSegmentIntersect(edges[static_cast<size_t>(l)], center, R))
1285  {
1286  toduplicate[static_cast<size_t>(l)].push_back(f[static_cast<size_t>(cur_facet)].vertices[static_cast<size_t>(i)]);
1287  added = true;
1288  }
1289  }
1290  }
1291  checked[static_cast<size_t>(f[static_cast<size_t>(cur_facet)].vertices[static_cast<size_t>(i)])] = true;
1292  if (added)
1293  {
1294  for (size_t j = 0; j < neigh.size(); ++j)
1295  tocheck.push(neigh[static_cast<size_t>(j)]);
1296  }
1297  }
1298  }
1299 }
1300 
1301 #ifdef RICH_MPI
1302 int Delaunay::findSomeOuterPoint(void)
1303 {
1304  const size_t cur_facet = Walk(olength);
1305  for (size_t i = 0; i < 3; ++i) {
1306  const int candidate = f.at(cur_facet).vertices[i];
1307  if (candidate < static_cast<int>(olength))
1308  return candidate;
1309  }
1310 #ifdef RICH_MPI
1311  int rank = 0;
1312  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1313  std::cout << "Problem with rank " << rank <<" chcked facet "<<cur_facet<<" olength "
1314  <<olength<< std::endl;
1315  throw UniversalError("Bad findsomeouterpoint");
1316 #endif
1317  assert(false && "something went wrong");
1318 }
1319 
1320 namespace
1321 {
1322  vector<int>
1323  calc_neighbors_own_edges
1324  (const Tessellation& t_proc,
1325  const vector<Edge>& edge_list, bool periodic)
1326  {
1327  vector<int> res;
1328  int rank;
1329  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1330  BOOST_FOREACH(const Edge& edge, edge_list)
1331  {
1332  const int other = (edge.neighbors.first + edge.neighbors.second) - rank;
1333  if (other < t_proc.GetPointNo())
1334  res.push_back(other);
1335  else
1336  if (periodic)
1337  res.push_back(t_proc.GetOriginalIndex(other));
1338  }
1339  if (periodic)
1340  {
1341  std::sort(res.begin(), res.end());
1342  res = unique(res);
1343  RemoveVal(res, rank);
1344  }
1345  return res;
1346  }
1347 
1348  stack<int> initialise_tocheck
1349  (const vector<int>& neightemp)
1350  {
1351  stack<int> res;
1352  for (size_t i = 0; i < neightemp.size(); ++i)
1353  res.push(neightemp[i]);
1354  return res;
1355  }
1356 
1357  vector<int> calc_self_intersection
1358  (const vector<Edge>& edge_list,
1359  const Circle& circle)
1360  {
1361  vector<int> res;
1362  for (size_t i = 0; i < edge_list.size(); ++i) {
1363  const Edge& edge = edge_list.at(i);
1364  if (edge_circle_intersect(edge, circle))
1365  res.push_back(static_cast<int>(i));
1366  }
1367  return res;
1368  }
1369 }
1370 
1371 vector<vector<int> > Delaunay::AddOuterFacetsMPI
1372 (int point,
1373  vector<vector<int> > &toduplicate,
1374  vector<int> &neigh,
1375  vector<bool> &checked,
1376  Tessellation const &tproc,
1377  const vector<Edge>& own_edges,
1378  bool periodic, std::vector<Vector2D> &periodic_add_self,
1379  std::vector<std::vector<Vector2D> >& periodic_add_others, Vector2D const &ll, Vector2D const &ur, bool recursive)
1380 {
1381  int rank;
1382  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1383  double procsize = tproc.GetWidth(rank);
1384  vector<vector<int> > res;
1385  vector<int> vtemp;
1386  std::vector<Vector2D> toadd;
1387  if (!recursive)
1388  {
1389  if (!periodic)
1390  res.resize(own_edges.size());
1391  else
1392  res.resize(1);
1393  }
1394  else
1395  if (periodic)
1396  res.resize(1);
1397  stack<int> tocheck = initialise_tocheck(FindContainingTetras(static_cast<int>(Walk(static_cast<size_t>(point))),
1398  point));
1399  if (recursive)
1400  {
1401  vector<int> allouter;
1402  for (size_t i = 0; i < toduplicate.size(); ++i)
1403  {
1404  for (size_t j = 0; j < toduplicate[i].size(); ++j)
1405  {
1406  vector<int> temp = FindContainingTetras
1407  (static_cast<int>(Walk(static_cast<size_t>(toduplicate[i][j]))), toduplicate[i][j]);
1408  for (size_t k = 0; k < temp.size(); ++k)
1409  allouter.push_back(temp[k]);
1410  }
1411  }
1412  sort(allouter.begin(), allouter.end());
1413  allouter = unique(allouter);
1414  for (size_t i = 0; i < allouter.size(); ++i)
1415  tocheck.push(allouter[i]);
1416  }
1417  while (!tocheck.empty())
1418  {
1419  int cur_facet = tocheck.top();
1420  tocheck.pop();
1421  for (size_t i = 0; i < 3; ++i)
1422  {
1423  bool added = false;
1424  int max_neigh = 0;
1425  if (f[static_cast<size_t>(cur_facet)].vertices[i] >= static_cast<int>(olength))
1426  continue;
1427  if (checked[static_cast<size_t>(f[static_cast<size_t>(cur_facet)].vertices[i])])
1428  continue;
1429  vector<int> neighs = FindContainingTetras(cur_facet, f[static_cast<size_t>(cur_facet)].vertices[i]);
1430 // int checking = f[static_cast<size_t>(cur_facet)].vertices[i];
1431  double R = 0;
1432  Vector2D center;
1433  for (size_t k = 0; k < neighs.size(); ++k)
1434  {
1435  center=GetCircleCenter(neighs[k], R);
1436  Circle circ(center,R);
1437  vector<int> cputosendto;
1438  if (recursive)
1439  find_affected_cells_recursive(tproc, rank, circ, cputosendto, toadd, periodic,ll,ur);
1440  else
1441  cputosendto = find_affected_cells(tproc, rank, circ, vtemp, periodic, toadd);
1442  if (!periodic)
1443  {
1444  // Remove self
1445  sort(cputosendto.begin(), cputosendto.end());
1446  cputosendto = unique(cputosendto);
1447  RemoveVal(cputosendto, rank);
1448  }
1449  if (periodic && recursive) // Remove self
1450  {
1451  if (!cputosendto.empty())
1452  {
1453  std::vector<int> toremove2;
1454  for (size_t l = 0; l < cputosendto.size(); ++l)
1455  if (cputosendto[l] == rank && fastabs(toadd[l]) < 0.01*procsize)
1456  toremove2.push_back(static_cast<int>(l));
1457  RemoveVector(cputosendto, toremove2);
1458  RemoveVector(toadd, toremove2);
1459  }
1460  }
1461  if (!recursive)
1462  {
1463  if (!periodic)
1464  {
1465  const vector<int> self_intersection = calc_self_intersection(own_edges, circ);
1466  if (!self_intersection.empty())
1467  added = true;
1468  BOOST_FOREACH(int sindex, self_intersection)
1469  res[static_cast<size_t>(sindex)].push_back(f[static_cast<size_t>(cur_facet)].vertices[i]);
1470  }
1471  }
1472  else
1473  {
1474  for (size_t jj = 0; jj < 3; ++jj)
1475  max_neigh = max(max_neigh, f[static_cast<size_t>(neighs[k])].vertices[jj]);
1476  }
1477  if (!cputosendto.empty())
1478  {
1479  added = true;
1480  for (size_t j = 0; j < cputosendto.size(); ++j)
1481  {
1482  if (cputosendto[j] == rank)
1483  {
1484  res[0].push_back(f[static_cast<size_t>(cur_facet)].vertices[i]);
1485  periodic_add_self.push_back(toadd[j]);
1486  continue;
1487  }
1488  size_t index = static_cast<size_t>(std::find(neigh.begin(), neigh.end(), cputosendto[j])
1489  - neigh.begin());
1490  if (index < neigh.size())
1491  {
1492  toduplicate.at(index).push_back(f[static_cast<size_t>(cur_facet)].vertices[i]);
1493  if (periodic)
1494  periodic_add_others.at(index).push_back(toadd[j]);
1495  }
1496  else
1497  {
1498  neigh.push_back(cputosendto[j]);
1499  toduplicate.push_back(vector<int>(1, f[static_cast<size_t>(cur_facet)].vertices[i]));
1500  if (periodic)
1501  periodic_add_others.push_back(vector<Vector2D>(1, toadd[j]));
1502  }
1503  }
1504  }
1505  }
1506  checked[static_cast<size_t>(f[static_cast<size_t>(cur_facet)].vertices[i])] = true;
1507  if (added || (recursive&&static_cast<size_t>(max_neigh) >= olength))
1508  {
1509  for (size_t j = 0; j < neighs.size(); ++j)
1510  tocheck.push(neighs[j]);
1511  }
1512  }
1513  }
1514  return res;
1515 }
1516 
1517 pair<vector<vector<int> >, vector<vector<int> > > Delaunay::findOuterPoints(const Tessellation& t_proc,
1518  const vector<Edge>& edge_list, const vector<Edge>& box_edges, vector<vector<int> > &NghostIndex, bool periodic,
1519  std::vector<int> &cpu_neigh, std::vector<Vector2D> &periodic_add_self, std::vector<std::vector<Vector2D> >
1520  &periodic_add_others, Vector2D const &ll, Vector2D const &ur)
1521 {
1522  cpu_neigh = calc_neighbors_own_edges(t_proc, edge_list, periodic);
1523  const size_t some_outer_point = findSomeOuterPoint();
1524  vector<vector<int> > to_duplicate(cpu_neigh.size());
1525  periodic_add_others.resize(to_duplicate.size());
1526  vector<bool> checked(olength, false);
1527  vector<vector<int> > self_points =
1528  AddOuterFacetsMPI
1529  (static_cast<int>(some_outer_point),
1530  to_duplicate, // indices of points to send
1531  cpu_neigh, // Rank of processes to send to
1532  checked,
1533  t_proc,
1534  box_edges, periodic, periodic_add_self, periodic_add_others,ll,ur);
1535  if (periodic)
1536  PeriodicGetRidDuplicatesSingle(self_points[0], periodic_add_self);
1537  for (size_t i = 0; i < to_duplicate.size(); ++i)
1538  {
1539  std::vector<size_t> sindex = sort_index(to_duplicate[i]);
1540  to_duplicate[i] = VectorValues(to_duplicate[i], sindex);
1541  if (periodic)
1542  {
1543  periodic_add_others[i] = VectorValues(periodic_add_others[i], sindex);
1544  PeriodicGetRidDuplicatesSingle(to_duplicate[i], periodic_add_others[i]);
1545  }
1546  else
1547  to_duplicate[i] = unique(to_duplicate[i]);
1548  }
1549  /*BOOST_FOREACH(vector<int>& line, to_duplicate)
1550  {
1551  sort(line.begin(), line.end());
1552  line = unique(line);
1553  }*/
1554 
1555  // Communication
1556  vector<vector<Vector2D> > incoming(cpu_neigh.size());
1557  vector<MPI_Request> req(cpu_neigh.size());
1558  vector<vector<double> > tosend(cpu_neigh.size());
1559  double dtemp = 0;
1560  for (size_t i = 0; i < cpu_neigh.size(); ++i)
1561  {
1562  const int dest = cpu_neigh[i];
1563  std::vector<Vector2D> cortosend = VectorValues(cor, to_duplicate[i]);
1564  if (periodic)
1565  {
1566  for (size_t j = 0; j < cortosend.size(); ++j)
1567  cortosend[j] += periodic_add_others[i][j];
1568  }
1569  tosend[i] = list_serialize(cortosend);
1570  int size = static_cast<int>(tosend[i].size());
1571  if (size == 0)
1572  MPI_Isend(&dtemp, 1, MPI_DOUBLE, dest, 1, MPI_COMM_WORLD, &req[i]);
1573  else
1574  {
1575  if (size < 2)
1576  throw UniversalError("Wrong send size");
1577  MPI_Isend(&tosend[i][0], size, MPI_DOUBLE, dest, 0, MPI_COMM_WORLD, &req[i]);
1578  }
1579  }
1580  for (size_t i = 0; i < cpu_neigh.size(); ++i)
1581  {
1582  vector<double> temprecv;
1583  MPI_Status status;
1584  MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
1585  int count;
1586  MPI_Get_count(&status, MPI_DOUBLE, &count);
1587  temprecv.resize(static_cast<size_t>(max(count, 1)));
1588  MPI_Recv(&temprecv[0], count, MPI_DOUBLE, status.MPI_SOURCE, status.MPI_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
1589  if (status.MPI_TAG == 0)
1590  {
1591  size_t location = static_cast<size_t>(std::find(cpu_neigh.begin(), cpu_neigh.end(),
1592  status.MPI_SOURCE) - cpu_neigh.begin());
1593  if (location >= cpu_neigh.size())
1594  throw UniversalError("Bad location in mpi exchange");
1595  try
1596  {
1597  incoming[location] = list_unserialize(temprecv, cor[0]);
1598  }
1599  catch (UniversalError &eo)
1600  {
1601  eo.AddEntry("Error in first send in triangulation", 0.0);
1602  throw;
1603  }
1604  }
1605  else
1606  if (status.MPI_TAG != 1)
1607  throw UniversalError("Wrong mpi tag");
1608  }
1609  if (req.size() > 0)
1610  MPI_Waitall(static_cast<int>(req.size()), &req[0], MPI_STATUSES_IGNORE);
1611  // Incorporate points recieved into triangulation
1612  if (!periodic)
1613  {
1614  for (size_t i = 0; i < self_points.size(); ++i)
1615  {
1616  std::vector<int> sindex;
1617  sort_index(self_points[i], sindex);
1618  self_points[i] = VectorValues(self_points[i], sindex);
1619  sindex = unique_index(self_points[i]);
1620  self_points[i] = VectorValues(self_points[i], sindex);
1621  }
1622  AddRigid(box_edges, self_points);
1623  }
1624  else
1625  {
1626  std::vector<size_t> indeces = sort_index(self_points[0]);
1627  self_points[0] = VectorValues(self_points[0], indeces);
1628  periodic_add_self = VectorValues(periodic_add_self, indeces);
1629  AddPeriodicMPI(self_points[0], periodic_add_self);
1630  }
1631  for (size_t i = 0; i < self_points.size(); ++i)
1632  for (size_t j = 0; j < self_points[i].size(); ++j)
1633  OrgIndex.push_back(self_points[i][j]);
1634 
1635  NghostIndex.clear();
1636  NghostIndex.resize(incoming.size());
1637  for (size_t i = 0; i < incoming.size(); ++i)
1638  {
1639  for (size_t j = 0; j < incoming.at(i).size(); ++j)
1640  {
1641  OrgIndex.push_back(static_cast<int>(cor.size() + j));
1642  NghostIndex[i].push_back(static_cast<int>(cor.size() + j));
1643  }
1644  AddBoundaryPoints(incoming.at(i));
1645  }
1646  MPI_Barrier(MPI_COMM_WORLD);
1647  return std::pair<vector<vector<int> >, vector<vector<int> > >
1648  (to_duplicate, self_points);
1649 }
1650 
1651 vector<vector<int> > Delaunay::boundary_intersection_check(const vector<Edge>& edges, const vector<vector<int> >& to_duplicate)
1652 {
1653  vector<vector<int> > res;
1654  res.reserve(edges.size());
1655  double R = 0;
1656  Vector2D center;
1657  BOOST_FOREACH(const Edge& edge, edges) {
1658  res.push_back(vector<int>());
1659  BOOST_FOREACH(const vector<int>& line, to_duplicate) {
1660  BOOST_FOREACH(const int index, line)
1661  {
1662  center=GetCircleCenter(index,R);
1663  const Circle circle
1664  (center, R);
1665  if (edge_circle_intersect(edge, circle))
1666  res.back().push_back(index);
1667  }
1668  }
1669  }
1670  return res;
1671 }
1672 
1673 pair<vector<vector<int> >, vector<int> > Delaunay::FindOuterPoints2
1674 (const Tessellation& t_proc,
1675  const vector<Edge>& edge_list,
1676  vector<vector<int> > &to_duplicate,
1677  vector<vector<int> >& self_points,
1678  const vector<Edge>& box_edges,
1679  vector<vector<int> > &NghostIndex,
1680  bool periodic,
1681  std::vector<int> &cpu_neigh, std::vector<Vector2D> &periodic_add_self, std::vector<std::vector<Vector2D> >
1682  &periodic_add_others, Vector2D const &ll, Vector2D const &ur)
1683 {
1684  vector<vector<int> > real_boundary_points(box_edges.size());
1685  if (!periodic)
1686  {
1687  const vector<vector<int> > boundary_points = boundary_intersection_check(box_edges, to_duplicate);
1688  for (size_t i = 0; i < boundary_points.size(); ++i)
1689  {
1690  sort(self_points.at(i).begin(), self_points.at(i).end());
1691  BOOST_FOREACH(int bp, boundary_points.at(i))
1692  {
1693  if (!binary_search(self_points.at(i).begin(),
1694  self_points.at(i).end(),
1695  bp))
1696  {
1697  real_boundary_points.at(i).push_back(bp);
1698  }
1699  }
1700  sort(real_boundary_points.at(i).begin(),
1701  real_boundary_points.at(i).end());
1702  real_boundary_points.at(i) = unique(real_boundary_points.at(i));
1703  }
1704  }
1705  else
1706  {
1707  for (size_t i = 0; i < self_points.size(); ++i)
1708  {
1709  std::vector<size_t> sindex = sort_index(self_points[i]);
1710  self_points[i] = VectorValues(self_points[i],sindex);
1711  periodic_add_self = VectorValues(periodic_add_self, sindex);
1712  }
1713  }
1714 
1715  vector<vector<int> > to_duplicate_2 = to_duplicate;
1716  vector<int> old_neighbors = cpu_neigh;
1717  //assert(!to_duplicate.empty());
1718  vector<bool> checked(olength, false);
1719  size_t some_outer_point = 0;
1720  if (!to_duplicate.empty())
1721  {
1722  if (!to_duplicate[0].empty())
1723  some_outer_point = to_duplicate[0][0];
1724  else
1725  some_outer_point = self_points.at(0).at(0);
1726  }
1727  else
1728  some_outer_point = self_points.at(0).at(0);
1729 
1730  std::vector<Vector2D> periodic_add_self2;
1731  std::vector<std::vector<Vector2D> > periodic_add_others2(periodic_add_others);
1732  std::vector<std::vector<int> > self_send = AddOuterFacetsMPI
1733  (static_cast<int>(some_outer_point),
1734  to_duplicate, // indices of points to send
1735  cpu_neigh, // Rank of processes to send to
1736  checked,
1737  t_proc,
1738  edge_list,
1739  periodic,
1740  periodic_add_self2,
1741  periodic_add_others2,ll,ur, true); // recursive
1742 
1743  // Communication
1744  int wsize;
1745  MPI_Comm_size(MPI_COMM_WORLD, &wsize);
1746  vector<int> totalk(static_cast<size_t>(wsize), 0);
1747  vector<int> scounts(totalk.size(), 1);
1748  for (size_t i = 0; i < cpu_neigh.size(); ++i)
1749  totalk[cpu_neigh[i]] = 1;
1750  int nrecv;
1751  MPI_Reduce_scatter(&totalk[0], &nrecv, &scounts[0], MPI_INT, MPI_SUM,
1752  MPI_COMM_WORLD);
1753 
1754  vector<MPI_Request> req(cpu_neigh.size());
1755  for (size_t i = 0; i < cpu_neigh.size(); ++i)
1756  MPI_Isend(&wsize, 1, MPI_INT, cpu_neigh[i], 3, MPI_COMM_WORLD, &req[i]);
1757  vector<int> talkwithme;
1758  for (int i = 0; i < nrecv; ++i)
1759  {
1760  MPI_Status status;
1761  MPI_Recv(&wsize, 1, MPI_INT, MPI_ANY_SOURCE, 3, MPI_COMM_WORLD, &status);
1762  talkwithme.push_back(status.MPI_SOURCE);
1763  }
1764  vector<size_t> indices;
1765  for (size_t i = 0; i < cpu_neigh.size(); ++i)
1766  {
1767  if (is_in(cpu_neigh[i], talkwithme))
1768  indices.push_back(i);
1769  }
1770 
1771  // Symmetrisation
1772  cpu_neigh = VectorValues(cpu_neigh, indices);
1773  to_duplicate = VectorValues(to_duplicate, indices);
1774  if (periodic)
1775  periodic_add_others2 = VectorValues(periodic_add_others2, indices);
1776  // Get rid of duplicate points
1777  int rank = 0;
1778  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1779  for (size_t i = 0; i < to_duplicate.size(); ++i)
1780  {
1781  if (periodic)
1782  {
1783  if (i < to_duplicate_2.size())
1784  PeriodicGetRidDuplicates(to_duplicate[i], periodic_add_others2[i], to_duplicate_2[i], periodic_add_others[i],t_proc.GetWidth(rank));
1785  else
1786  PeriodicGetRidDuplicatesSingle(to_duplicate[i], periodic_add_others2[i]);
1787  }
1788  else
1789  {
1790  sort(to_duplicate[i].begin(), to_duplicate[i].end());
1791  to_duplicate[i] = unique(to_duplicate[i]);
1792  }
1793  }
1794  vector<vector<int> > messages(to_duplicate.size());
1795  if (!periodic)
1796  {
1797  for (size_t i = 0; i < to_duplicate.size(); ++i)
1798  {
1799  const vector<int>::const_iterator it = find(old_neighbors.begin(), old_neighbors.end(), cpu_neigh.at(i));
1800  for (size_t j = 0; j < to_duplicate.at(i).size(); ++j)
1801  {
1802  if (it != old_neighbors.end())
1803  {
1804  const size_t my_index = static_cast<size_t>(it - old_neighbors.begin());
1805  if (!binary_search(to_duplicate_2.at(my_index).begin(), to_duplicate_2.at(my_index).end(),
1806  to_duplicate.at(i).at(j)))
1807  messages.at(i).push_back(to_duplicate.at(i).at(j));
1808  }
1809  else
1810  messages.at(i).push_back(to_duplicate.at(i).at(j));
1811  }
1812  }
1813  }
1814  if (req.size() > 0)
1815  MPI_Waitall(static_cast<int>(cpu_neigh.size()), &req[0], MPI_STATUSES_IGNORE);
1816  MPI_Barrier(MPI_COMM_WORLD);
1817  // Point exchange
1818  req.clear();
1819  req.resize(cpu_neigh.size());
1820  double dtemp = 0;
1821  vector<vector<Vector2D> > incoming(cpu_neigh.size());
1822  vector<vector<double> > tosend(cpu_neigh.size());
1823  for (size_t i = 0; i < cpu_neigh.size(); ++i)
1824  {
1825  const int dest = cpu_neigh[i];
1826  if (periodic)
1827  {
1828  std::vector<Vector2D> corsend = VectorValues(cor, to_duplicate[i]);
1829  for (size_t j = 0; j < corsend.size(); ++j)
1830  corsend[j] += periodic_add_others2[i][j];
1831  tosend[i] = list_serialize(corsend);
1832  }
1833  else
1834  tosend[i] = list_serialize(VectorValues(cor, messages[i]));
1835  int size = static_cast<int>(tosend[i].size());
1836  if (size > 0)
1837  {
1838  if (size < 2)
1839  throw UniversalError("Wrong send size");
1840  MPI_Isend(&tosend[i][0], size, MPI_DOUBLE, dest, 0, MPI_COMM_WORLD, &req[i]);
1841  }
1842  else
1843  MPI_Isend(&dtemp, 1, MPI_DOUBLE, dest, 1, MPI_COMM_WORLD, &req[i]);
1844  }
1845  for (size_t i = 0; i < cpu_neigh.size(); ++i)
1846  {
1847  vector<double> temprecv;
1848  MPI_Status status;
1849  MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
1850  int count;
1851  MPI_Get_count(&status, MPI_DOUBLE, &count);
1852  temprecv.resize(static_cast<size_t>(count));
1853  MPI_Recv(&temprecv[0], count, MPI_DOUBLE, status.MPI_SOURCE, status.MPI_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
1854  if (status.MPI_TAG == 0)
1855  {
1856  size_t location = static_cast<size_t>(std::find(cpu_neigh.begin(), cpu_neigh.end(), status.MPI_SOURCE) -
1857  cpu_neigh.begin());
1858  if (location >= cpu_neigh.size())
1859  throw UniversalError("Bad location in mpi exchange");
1860  try
1861  {
1862  incoming[location] = list_unserialize(temprecv, cor[0]);
1863  }
1864  catch (UniversalError &eo)
1865  {
1866  eo.AddEntry("Error in second send in triangulation", 0.0);
1867  eo.AddEntry("Mpi status", static_cast<double>(status.MPI_SOURCE));
1868  eo.AddEntry("Mpi tag", static_cast<double>(status.MPI_TAG));
1869  eo.AddEntry("Mpi count", static_cast<double>(count));
1870  throw;
1871  }
1872  }
1873  else
1874  if (status.MPI_TAG != 1)
1875  throw UniversalError("Wrong mpi tag");
1876  }
1877  if (req.size() > 0)
1878  MPI_Waitall(static_cast<int>(req.size()), &req[0], MPI_STATUSES_IGNORE);
1879  MPI_Barrier(MPI_COMM_WORLD);
1880  // Incorporate points recieved into triangulation
1881  NghostIndex.resize(incoming.size());
1882  for (size_t i = 0; i < incoming.size(); ++i)
1883  {
1884  for (size_t j = 0; j < incoming.at(i).size(); ++j)
1885  {
1886  NghostIndex[i].push_back(static_cast<int>(cor.size() + j));
1887  OrgIndex.push_back(static_cast<int>(cor.size() + j));
1888  }
1889  AddBoundaryPoints(incoming.at(i));
1890  }
1891 
1892  // ADD SELF SEND FROM SECOND RUN IN PERIODIC
1893  if (periodic)
1894  {
1895  PeriodicGetRidDuplicates(self_send[0], periodic_add_self2, self_points[0], periodic_add_self,t_proc.GetWidth(rank));
1896  AddPeriodicMPI(self_send[0], periodic_add_self2);
1897  for (size_t j = 0; j < self_send[0].size(); ++j)
1898  OrgIndex.push_back(self_send[0][j]);
1899  // Add the dulicates from first run
1900  for (size_t i = 0; i < to_duplicate_2.size(); ++i)
1901  to_duplicate[i].insert(to_duplicate[i].begin(), to_duplicate_2[i].begin(), to_duplicate_2[i].end());
1902  }
1903  else
1904  {
1905  AddRigid(box_edges, real_boundary_points);
1906  for (size_t i = 0; i < real_boundary_points.size(); ++i)
1907  for (size_t j = 0; j < real_boundary_points[i].size(); ++j)
1908  OrgIndex.push_back(real_boundary_points[i][j]);
1909 
1910  for (size_t i = 0; i < self_points.size(); ++i)
1911  if (!real_boundary_points[i].empty())
1912  self_points[i].insert(self_points[i].end(), real_boundary_points[i].begin(),
1913  real_boundary_points[i].end());
1914  }
1915  return std::pair<vector<vector<int> >, vector<int> >(to_duplicate, cpu_neigh);
1916 }
1917 
1918 std::pair<vector<vector<int> >, vector<int> > Delaunay::BuildBoundary
1919 (OuterBoundary const* obc,
1920  Tessellation const& tproc,
1921  vector<vector<int> >& Nghost)
1922 {
1923  vector<Edge> edges;
1924  OrgIndex.clear();
1925  int rank;
1926  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1927  vector<int> edge_index = tproc.GetCellEdges(rank);
1928  /*
1929  if(rank==0){
1930  ofstream dump("edge_index.txt");
1931  for(size_t i=0;i<edge_index.size();++i)
1932  dump << edge_index.at(i) << endl;
1933  dump.close();
1934  assert(false);
1935  }
1936  */
1937  for (size_t i = 0; i < edge_index.size(); ++i)
1938  edges.push_back(tproc.GetEdge(edge_index[i]));
1939  bool periodic = obc->GetBoundaryType() == Periodic;
1940  vector<Edge> box_edges;
1941  if (periodic)
1942  box_edges = edges;
1943  else
1944  box_edges = obc->GetBoxEdges();
1945  std::vector<int> cpu_neigh;
1946  std::vector<Vector2D> periodic_add_self;
1947  std::vector<std::vector<Vector2D> > periodic_add_others;
1948  std::pair<vector<vector<int> >, vector<vector<int> > > to_duplicate =
1949  findOuterPoints(tproc, edges, box_edges, Nghost, periodic, cpu_neigh, periodic_add_self, periodic_add_others,
1950  Vector2D(obc->GetGridBoundary(Left), obc->GetGridBoundary(Down)),
1951  Vector2D(obc->GetGridBoundary(Right), obc->GetGridBoundary(Up)));
1952  std::pair<vector<vector<int> >, vector<int> > res = FindOuterPoints2(tproc, edges, to_duplicate.first, to_duplicate.second, box_edges, Nghost, periodic, cpu_neigh,
1953  periodic_add_self, periodic_add_others, Vector2D(obc->GetGridBoundary(Left), obc->GetGridBoundary(Down)),
1954  Vector2D(obc->GetGridBoundary(Right), obc->GetGridBoundary(Up)));
1955  radius.resize(f.size());
1956  int n = int(f.size());
1957  for (int i = 0; i < n; ++i)
1958  radius[static_cast<size_t>(i)] = CalculateRadius(i);
1959  return res;
1960 }
1961 #endif // RICH_MPI
1962 
1963 int Delaunay::GetOrgIndex(int index)const
1964 {
1965  if (index < static_cast<int>(olength))
1966  return static_cast<int>(olength);
1967  else
1968  return OrgIndex.at(static_cast<size_t>(index) - 3 - olength);
1969 }
unique
vector< T > unique(vector< T > const &v)
Returns a vector containing only unique elements.
Definition: utils.hpp:376
Delaunay::get_cor
double get_cor(int index, int dim) const
Returns a coordinate.
Definition: Delaunay.cpp:709
TripleConstRef::third
const T & third
Reference to third item.
Definition: triplet.hpp:23
VectorValues
vector< T > VectorValues(vector< T > const &v, vector< int > const &index)
Returns only the values with indeces in index.
Definition: utils.hpp:326
Delaunay::Changelength
void Changelength(int n)
Changes the cor length.
Definition: Delaunay.cpp:676
incircle
double incircle(const Vector2D &point_1, const Vector2D &point_2, const Vector2D &point_3, const Vector2D &point_4)
Checks whether the 4th point is inside, outside or on the counterclockwise circle created by the firs...
Definition: geotests.cpp:567
Delaunay::AddBoundaryPoints
void AddBoundaryPoints(vector< Vector2D > const &points)
Adds the points to the tessellation, used for boundary points.
Definition: Delaunay.cpp:754
Delaunay::AddAditionalPoint
void AddAditionalPoint(Vector2D const &vec)
Adds a point to the cor vector. Used in periodic boundaries with AMR.
Definition: Delaunay.cpp:767
sort_index
void sort_index(const vector< T > &arr, vector< int > &res)
Returns the indeces of a sort.
Definition: utils.hpp:500
Parallel
Vector2D Parallel(Edge const &edge)
Calculates a unit vector parallel to an edge.
Definition: Edge.cpp:83
Delaunay::set_point
void set_point(int index, Vector2D p)
Change Mesh point.
Definition: Delaunay.cpp:734
TripleConstRef::second
const T & second
Reference to second item.
Definition: triplet.hpp:20
RemoveVector
void RemoveVector(vector< T > &v, vector< int > &indeces)
Removes the elements in v given by indeces.
Definition: utils.hpp:275
Tessellation
Abstract class for tessellation.
Definition: tessellation.hpp:21
Tessellation::GetPointNo
virtual int GetPointNo(void) const =0
Get Total number of mesh generating points.
Delaunay::get_facet
const facet & get_facet(int index) const
Returns a facet.
Definition: Delaunay.cpp:691
Tessellation::GetOriginalIndex
virtual int GetOriginalIndex(int point) const
Returns the original index of the duplicated point.
Definition: tessellation.cpp:5
facet
Delanauy triangle data structure. Keeps the indexes of the vertices (right hand fashion) and the neig...
Definition: facet.hpp:19
UniversalError
Container for error reports.
Definition: universal_error.hpp:14
Delaunay::update
void update(const vector< Vector2D > &points, std::vector< std::pair< Vector2D, Vector2D > >const &cpoints)
Updates the triangulation.
Definition: Delaunay.cpp:504
Delaunay::getCor
const vector< Vector2D > & getCor(void) const
Access to coordinates.
Definition: Delaunay.cpp:686
facet::vertices
Triplet< int > vertices
Indices of vertices.
Definition: facet.hpp:24
Delaunay::get_length
int get_length(void) const
Returns the number of points.
Definition: Delaunay.cpp:724
Triplet::set
void set(const T &first_i, const T &second_i, const T &third_i)
Changle all items.
Definition: triplet.hpp:77
Delaunay::GetTotalLength
int GetTotalLength(void)
Returns the length of all the points (included duplicated)
Definition: Delaunay.cpp:749
Edge
Interface between two cells.
Definition: Edge.hpp:13
Delaunay::GetCircleCenter
Vector2D GetCircleCenter(int index, double &R) const
Returns the center of the circumscribed circle of a facet.
Definition: Delaunay.cpp:1202
max
double max(vector< double > const &v)
returns the maximal term in a vector
Definition: utils.cpp:52
Vector2D::y
double y
Component in the y direction.
Definition: geometry.hpp:48
TripleConstRef
A collection of three identical references.
Definition: triplet.hpp:12
Delaunay::GetMeshPoints
vector< Vector2D > & GetMeshPoints(void)
Returns a refrence to the points.
Definition: Delaunay.cpp:744
Tessellation::GetEdge
virtual Edge const & GetEdge(int index) const =0
Returns edge (interface between cells)
Delaunay::BuildBoundary
void BuildBoundary(OuterBoundary const *obc, vector< Edge > const &edges)
Builds the boundary points.
Definition: Delaunay.cpp:1167
Delaunay::~Delaunay
~Delaunay(void)
Default destructor.
Definition: Delaunay.cpp:177
Delaunay::get_num_facet
int get_num_facet(void) const
Returns the number of facets.
Definition: Delaunay.cpp:719
Delaunay
The Delaunay data structure. Gets a set of points and constructs the Delaunay tessellation.
Definition: Delaunay.hpp:30
ScalarProd
double ScalarProd(Vector3D const &v1, Vector3D const &v2)
Scalar product of two vectors.
Definition: Vector3D.cpp:185
RemoveVal
void RemoveVal(vector< T > &vec, T val)
Removes the first occurence of val inside a vector.
Definition: utils.hpp:431
OuterBoundary::GetBoundaryType
virtual BoundaryType GetBoundaryType(void) const =0
Returns the boundary type.
Delaunay::GetCorSize
int GetCorSize(void) const
Gets the size of the cor vector.
Definition: Delaunay.cpp:772
delaunay_loggers::DelaunayLogger::output
virtual void output(vector< Vector2D > const &cor, vector< facet > const &f)
Dumps output.
Definition: delaunay_logger.cpp:9
Circle
A circle.
Definition: shape_2d.hpp:27
Edge::vertices
std::pair< Vector2D, Vector2D > vertices
Points at the ends of the edge.
Definition: Edge.hpp:18
ReArrangeVector
void ReArrangeVector(vector< T > &v, vector< int > const &indeces)
Rearranges the vector according to the indeces.
Definition: utils.hpp:477
Tessellation::GetWidth
virtual double GetWidth(int index) const =0
Returns the effective width of a cell.
Delaunay::GetOrgIndex
int GetOrgIndex(int index) const
Retrieves the original index of a point (in case a point was duplicated)
Definition: Delaunay.cpp:1963
Delaunay::get_point
Vector2D get_point(size_t index) const
Returns a point.
Definition: Delaunay.cpp:704
find_affected_cells
vector< int > find_affected_cells(const Tessellation &tess, int index, Circle &circle, vector< int > &vtemp, bool periodic, std::vector< Vector2D > &periodic_add)
Non recursive version of find affected cells. Only searches one degree of separation.
Definition: find_affected_cells.cpp:136
min
double min(vector< double > const &v)
Returns the minimal term in a vector.
Definition: utils.cpp:44
Delaunay::logger
delaunay_loggers::DelaunayLogger * logger
Diagnostics.
Definition: Delaunay.hpp:277
OuterBoundary::GetBoxEdges
vector< Edge > GetBoxEdges(void) const
Returns the outer box as a set of edges in the order: Right, Up, Left and Down. All neighbors are set...
Definition: OuterBoundary.cpp:5
HilbertOrder
vector< int > HilbertOrder(vector< Vector2D > const &cor, int num, int innernum=0)
Returns the Hilber curve ordering.
Definition: HilbertOrder.cpp:103
facet::neighbors
Triplet< int > neighbors
Indices of neighboring facets.
Definition: facet.hpp:27
Vector2D::Set
void Set(double ix, double iy)
Set vector components.
Definition: geometry.cpp:39
find_affected_cells_recursive
void find_affected_cells_recursive(const Tessellation &tess, int index, const Circle &circle, vector< int > &res, std::vector< Vector2D > &added, bool periodic, Vector2D const &ll, Vector2D const &ur)
Recursively finds all cells that intersect a circle.
Definition: find_affected_cells.cpp:174
unique_index
vector< int > unique_index(vector< T > const &v)
Returns a vector containing only indeces of unique elements.
Definition: utils.hpp:397
Delaunay::get_facet_coordinate
double get_facet_coordinate(int Facet, int vertice, int dim)
Returns a coordinate of a vertice.
Definition: Delaunay.cpp:696
Edge::neighbors
std::pair< int, int > neighbors
Neighboring cells.
Definition: Edge.hpp:21
Triplet::third
T third
Third item.
Definition: triplet.hpp:50
Delaunay::ChangeOlength
void ChangeOlength(int n)
Changes the cor olength.
Definition: Delaunay.cpp:671
edge_circle_intersect
bool edge_circle_intersect(const Edge &edge, const Circle &circle)
Determines if an edge and a circle intersect.
Definition: find_affected_cells.cpp:63
TripleConstRef::first
const T & first
Reference to first item.
Definition: triplet.hpp:17
Delaunay::GetOriginalLength
int GetOriginalLength(void) const
Returns the original length of the points (without duplicated points)
Definition: Delaunay.cpp:739
abs
double abs(Vector3D const &v)
Norm of a vector.
Definition: Vector3D.cpp:44
UniversalError::AddEntry
void AddEntry(std::string const &field, double value)
Adds an entry to the list.
Definition: universal_error.cpp:20
Tessellation::GetCellEdges
virtual vector< int > const & GetCellEdges(int index) const =0
Returns the indexes of a cell&#39;s edges.
Delaunay::get_last_loc
int get_last_loc(void) const
Returns the last location, a number used to identify the fact that the neighbor of a facet is empty...
Definition: Delaunay.cpp:729
OuterBoundary
Abstract class for geometric boundary conditions for the tessellation.
Definition: OuterBoundary.hpp:19
Delaunay::triangle_area
double triangle_area(int index)
Returns the area of the triangle. Negative result means the triangle isn&#39;t right handed.
Definition: Delaunay.cpp:491
Triplet::second
T second
Second item.
Definition: triplet.hpp:47
Delaunay.hpp
Delaunay triangulation with mpi.
Vector2D
2D Mathematical vector
Definition: geometry.hpp:15
Delaunay::build_delaunay
void build_delaunay(vector< Vector2D >const &vp, std::vector< std::pair< Vector2D, Vector2D > >const &cpoints)
Builds the Delaunay tessellation.
Definition: Delaunay.cpp:376
Delaunay::GetFacetRadius
double GetFacetRadius(int facet) const
return the facet&#39;s radius
Definition: Delaunay.cpp:666
OuterBoundary::GetGridBoundary
virtual double GetGridBoundary(Directions dir) const =0
Returns the boundary coordinate.
orient2d
double orient2d(const TripleConstRef< Vector2D > &points)
Checks whether 3 given points are on a counterclockwise circle, clockwise circle or colinear...
Definition: geotests.cpp:76
fastabs
double fastabs(Vector2D const &v)
Norm of a vector, less accurate.
Definition: geometry.cpp:24
Vector2D::x
double x
Component in the x direction.
Definition: geometry.hpp:45
Delaunay::Delaunay
Delaunay(void)
Class constructor.
Definition: Delaunay.cpp:148
Delaunay::ChangeCor
vector< Vector2D > & ChangeCor(void)
Allows to change the cor.
Definition: Delaunay.cpp:681
Triplet::first
T first
First item.
Definition: triplet.hpp:44
Delaunay::GetOriginalIndex
int GetOriginalIndex(int NewPoint) const
Returns the original index of the duplicated point in Periodic Boundary conditions.
Definition: Delaunay.cpp:661
Generated on Thu Jun 4 2020 18:12:17 by   doxygen 1.8.13