VoronoiMesh.cpp
1 #include "VoronoiMesh.hpp"
2 #include <cmath>
3 #include "../misc/simple_io.hpp"
4 #include "hdf5_logger.hpp"
5 #ifdef RICH_MPI
6 #include <mpi.h>
7 #include "../mpi/mpi_commands.hpp"
8 #endif
9 
10 using std::abs;
11 using std::endl;
12 using std::cout;
13 
14 namespace
15 {
16 #ifdef RICH_MPI
17  template<class T> void tidy(vector<T>& v)
18  {
19  if (!v.empty()) {
20  sort(v.begin(), v.end());
21  v = unique(v);
22  }
23  }
24 #endif
25 
26  vector<Vector2D> UpdatePoints(vector<Vector2D> const& points, OuterBoundary const* obc)
27  {
28  if (obc->GetBoundaryType() == Rectengular)
29  return points;
30  vector<Vector2D> res;
31  res.reserve(points.size());
32  int npoints = static_cast<int>(points.size());
33  const double dx = obc->GetGridBoundary(Right) - obc->GetGridBoundary(Left);
34  const double dy = obc->GetGridBoundary(Up) - obc->GetGridBoundary(Down);
35  for (int i = 0; i < npoints; ++i)
36  {
37  Vector2D temp(points[static_cast<size_t>(i)]);
38  if (obc->GetBoundaryType() == Periodic)
39  {
40  if (temp.x > obc->GetGridBoundary(Right))
41  temp.x -= dx;
42  if (temp.x < obc->GetGridBoundary(Left))
43  temp.x += dx;
44  if (temp.y > obc->GetGridBoundary(Up))
45  temp.y -= dy;
46  if (temp.y < obc->GetGridBoundary(Down))
47  temp.y += dy;
48  }
49  if (obc->GetBoundaryType() == HalfPeriodic)
50  {
51  if (temp.x > obc->GetGridBoundary(Right))
52  temp.x -= dx;
53  if (temp.x < obc->GetGridBoundary(Left))
54  temp.x += dx;
55  }
56  res.push_back(temp);
57  }
58  return res;
59  }
60 
61  void CheckInput(std::vector<std::pair<Vector2D, Vector2D> > convexhull, std::vector<Vector2D> const& cor)
62  {
63  size_t N = cor.size();
64  for (size_t i = 0; i < N; ++i)
65  {
66  if (!PointInCell(convexhull, cor[i]))
67  {
68  std::cout << "Point not in cell in checkinput" << std::endl;
69  std::cout << "Bad point " << cor[i].x << ", " << cor[i].y << std::endl;
70  std::cout << "Bounding cell " << std::endl;
71  for (size_t j = 0; j < convexhull.size(); ++j)
72  {
73  std::cout << convexhull[j].first.x << ", " << convexhull[j].first.y <<" "<<convexhull[j].second.x << ", " <<
74  convexhull[j].second.y << std::endl;
75  }
76  UniversalError eo("Bad check point");
77  eo.AddEntry("Point index", static_cast<double>(i));
78  throw eo;
79  }
80  }
81  }
82 }
83 
84 VoronoiMesh::VoronoiMesh
85 (vector<Vector2D> const& points,
86  OuterBoundary const& bc, bool HOrder) :
87  logger(0),
88  eps(1e-8),
89  obc(0),
90  cell_edges(vector<Edge>()),
91  edges(vector<Edge>()),
92  CM(vector<Vector2D>()),
93  mesh_vertices(vector<vector<int> >()),
94  Tri(),
95  GhostProcs(vector<int>()),
96  GhostPoints(vector<vector<int> >()),
97  SentProcs(vector<int>()),
98  SentPoints(vector<vector<int> >()),
99  selfindex(vector<size_t>()),
100  NGhostReceived(vector<vector<int> >()),
101  OrgCorner(),
102  Nextra(0)
103 {
104  Initialise(points, &bc, HOrder);
105 }
106 
107 #ifdef RICH_MPI
108 VoronoiMesh::VoronoiMesh
109 (Tessellation const& proctess,
110  vector<Vector2D> const& points,
111  OuterBoundary const& bc,
112  bool HOrder) :
113  logger(0),
114  eps(1e-8),
115  obc(0),
116  cell_edges(vector<Edge>()),
117  edges(vector<Edge>()),
118  CM(vector<Vector2D>()),
119  mesh_vertices(vector<vector<int> >()),
120  Tri(),
121  GhostProcs(vector<int>()),
122  GhostPoints(vector<vector<int> >()),
123  SentProcs(vector<int>()),
124  SentPoints(vector<vector<int> >()),
125  selfindex(vector<size_t>()),
126  NGhostReceived(vector<vector<int> >()),
127  OrgCorner(),
128  Nextra(0)
129 {
130  Initialise(points, proctess, &bc, HOrder);
131 }
132 #endif
133 
134 vector<int> VoronoiMesh::AddPointsAlongEdge(size_t point, vector<vector<int> > const&copied,
135  int side)
136 {
137  int ncopy = static_cast<int>(copied[static_cast<size_t>(side)].size());
138  Vector2D vec = Tri.get_point(point);
139  vector<double> dist(static_cast<size_t>(ncopy));
140  for (size_t i = 0; i < copied[static_cast<size_t>(side)].size(); ++i)
141  dist[i] = vec.distance(Tri.get_point(static_cast<size_t>(copied[static_cast<size_t>(side)][i])));
142  const int copylength = min(7, static_cast<int>(copied[static_cast<size_t>(side)].size()) - 1);
143  vector<int> index, toadd(static_cast<size_t>(copylength));
144  sort_index(dist, index);
145  for (int i = 0; i < copylength; ++i)
146  toadd[static_cast<size_t>(i)] = copied[static_cast<size_t>(side)][static_cast<size_t>(index[static_cast<size_t>(i) + 1])];
147  return toadd;
148 }
149 
150 Vector2D VoronoiMesh::CalcFaceVelocity(Vector2D wl, Vector2D wr, Vector2D rL, Vector2D rR,
151  Vector2D f)const
152 {
153  const Vector2D wprime = ScalarProd(wl - wr, f - (rR + rL) / 2)*(rR - rL) / pow(abs(rR - rL), 2);
154  return 0.5*(wl + wr) + wprime;
155 }
156 
157 bool VoronoiMesh::NearBoundary(int index) const
158 {
159  const int n = int(mesh_vertices[static_cast<size_t>(index)].size());
160  const int N = Tri.get_length();
161  for (int i = 0; i < n; ++i)
162  {
163  const int n0 = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][static_cast<size_t>(i)])].neighbors.first;
164  const int n1 = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][static_cast<size_t>(i)])].neighbors.second;
165  if (n0 < 0 || n1 < 0 || n0 >= N || n1 >= N)
166  return true;
167  }
168  return false;
169 }
170 
171 VoronoiMesh::~VoronoiMesh(void) {}
172 
173 int VoronoiMesh::GetOriginalIndex(int point) const
174 {
175  int npoints = GetPointNo();
176  if (point < npoints)
177  return point;
178  else
179  return Tri.GetOrgIndex(point);
180 }
181 
182 vector<size_t> VoronoiMesh::GetSelfPoint(void)const
183 {
184  return selfindex;
185 }
186 
187 VoronoiMesh::VoronoiMesh(void) :
188  logger(0),
189  eps(1e-8),
190  obc(0),
191  cell_edges(vector<Edge>()),
192  edges(vector<Edge>()),
193  CM(vector<Vector2D>()),
194  mesh_vertices(vector<vector<int> >()),
195  Tri(),
196  GhostProcs(vector<int>()),
197  GhostPoints(vector<vector<int> >()),
198  SentProcs(vector<int>()),
199  SentPoints(vector<vector<int> >()),
200  selfindex(vector<size_t>()),
201  NGhostReceived(vector<vector<int> >()),
202  OrgCorner(),
203  Nextra(0) {}
204 
205 VoronoiMesh::VoronoiMesh(VoronoiMesh const& other) :
206  logger(other.logger),
207  eps(other.eps),
208  obc(other.obc),
209  cell_edges(other.cell_edges),
210  edges(other.edges),
211  CM(other.CM),
212  mesh_vertices(other.mesh_vertices),
213  Tri(other.Tri),
214  GhostProcs(other.GhostProcs),
215  GhostPoints(other.GhostPoints),
216  SentProcs(other.SentProcs),
217  SentPoints(other.SentPoints),
218  selfindex(other.selfindex),
219  NGhostReceived(other.NGhostReceived),
220  OrgCorner(),
221  Nextra(other.Nextra)
222 {}
223 
224 void VoronoiMesh::build_v()
225 {
226  Vector2D center, center_temp;
227  int j;
228  facet to_check;
229  Edge edge_temp;
230  Vector2D p_temp;
231  mesh_vertices.clear();
232  mesh_vertices.resize(static_cast<size_t>(Tri.get_length()));
233  edges.reserve(static_cast<size_t>(Tri.get_length()*3.5));
234  int N = Tri.GetOriginalLength();
235  for (int i = 0; i < N; ++i)
236  mesh_vertices[static_cast<size_t>(i)].reserve(7);
237  int Nfacets = Tri.get_num_facet();
238  vector<Vector2D> centers(static_cast<size_t>(Nfacets));
239  double R = 0;
240  for (int i = 0; i < Nfacets; ++i)
241  centers[static_cast<size_t>(i)] = Tri.GetCircleCenter(i,R);
242  for (int i = 0; i < Nfacets; ++i)
243  {
244  center = centers[static_cast<size_t>(i)];
245  to_check = Tri.get_facet(i);
246  for (j = 0; j < 3; ++j)
247  {
248  if (to_check.neighbors[static_cast<size_t>(j)] == Tri.get_last_loc())
249  continue;
250  if (to_check.neighbors[static_cast<size_t>(j)] < i)
251  continue;
252  center_temp = centers[static_cast<size_t>(to_check.neighbors[static_cast<size_t>(j)])];
253  {
254  edge_temp.vertices.first = center;
255  edge_temp.vertices.second = center_temp;
256  edge_temp.neighbors.first = to_check.vertices[static_cast<size_t>(j)];
257  edge_temp.neighbors.second = to_check.vertices[static_cast<size_t>(j + 1) % 3];
258  if (legal_edge(&edge_temp))
259  {
260  // I added a change here, if edge has zero length I don't add it.
261  if (edge_temp.GetLength() > eps*sqrt(Tri.GetFacetRadius(i)*
262  Tri.GetFacetRadius(to_check.neighbors[static_cast<size_t>(j)])))
263  {
264  {
265  if (edge_temp.neighbors.first < Tri.GetOriginalLength())
266  mesh_vertices[static_cast<size_t>(edge_temp.neighbors.first)].push_back(static_cast<int>(edges.size()));
267  if (edge_temp.neighbors.second < Tri.GetOriginalLength())
268  mesh_vertices[static_cast<size_t>(edge_temp.neighbors.second)].push_back(static_cast<int>(edges.size()));
269  edges.push_back(edge_temp);
270  }
271  }
272  }
273  }
274  }
275  }
276 }
277 
278 namespace
279 {
280  std::vector<std::pair<Vector2D, Vector2D> > calc_procpoints(const OuterBoundary& bc)
281  {
282  std::vector<std::pair<Vector2D, Vector2D> > res(4);
283  res[0] = std::pair<Vector2D, Vector2D>(Vector2D(bc.GetGridBoundary(Left), bc.GetGridBoundary(Down)),
284  Vector2D(bc.GetGridBoundary(Right), bc.GetGridBoundary(Down)));
285  res[1] = std::pair<Vector2D, Vector2D>(Vector2D(bc.GetGridBoundary(Right), bc.GetGridBoundary(Down)),
286  Vector2D(bc.GetGridBoundary(Right), bc.GetGridBoundary(Up)));
287  res[2] = std::pair<Vector2D, Vector2D>(Vector2D(bc.GetGridBoundary(Right), bc.GetGridBoundary(Up)),
288  Vector2D(bc.GetGridBoundary(Left), bc.GetGridBoundary(Up)));
289  res[2] = std::pair<Vector2D, Vector2D>(Vector2D(bc.GetGridBoundary(Left), bc.GetGridBoundary(Up)),
290  Vector2D(bc.GetGridBoundary(Left), bc.GetGridBoundary(Down)));
291  return res;
292  }
293 
294  /*
295  Vector2D GetReflection(OuterBoundary const& bc, size_t index, Vector2D const& point)
296  {
297  switch (index)
298  {
299  case 0:
300  return point + Vector2D(2 * (bc.GetGridBoundary(Right) - point.x), 0);
301  case 1:
302  return point + Vector2D(0, 2 * (bc.GetGridBoundary(Up) - point.y));
303  case 2:
304  return point + Vector2D(2 * (bc.GetGridBoundary(Left) - point.x), 0);
305  case 3:
306  return point + Vector2D(0, 2 * (bc.GetGridBoundary(Down) - point.y));
307  }
308  throw UniversalError("Wrong index in VoronoiMesh::GetReflection");
309  }*/
310 }
311 
312 void VoronoiMesh::Initialise(vector<Vector2D>const& pv, OuterBoundary const* _bc, bool reorder)
313 {
314  obc = _bc;
315  vector<Vector2D> points;
316  if (reorder)
317  points = VectorValues(pv, HilbertOrder(pv, static_cast<int>(pv.size())));
318  else
319  points = pv;
320  // Check that points are all isnide domain
321  std::vector<std::pair<Vector2D, Vector2D> > cpoints = calc_procpoints(*obc);
322  points = UpdatePoints(points, obc);
323  CheckInput(cpoints, points);
324 
325  Tri.build_delaunay(points, cpoints);
326 
327  Nextra = static_cast<int>(Tri.ChangeCor().size());
328  Tri.BuildBoundary(_bc, _bc->GetBoxEdges());
329 
330  eps = 1e-8;
331  edges.clear();
332  GhostPoints.clear();
333  GhostProcs.clear();
334  build_v();
335 
336  if (logger)
337  logger->output(*this);
338 
339  CM.resize(Tri.getCor().size());
340  for (size_t i = 0; i < pv.size(); ++i)
341  CM[i] = CalcCellCM(i);
342 
343  size_t counter = pv.size() + 3;
344  if (_bc->GetBoundaryType() == Periodic)
345  {
346  for (int i = static_cast<int>(counter); i < Tri.GetCorSize(); ++i)
347  {
348  int NorgIndex = Tri.GetOrgIndex(i);
349  if (NorgIndex < Nextra)
350  {
351  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] +
352  (Tri.get_point(static_cast<size_t>(i)) -
353  Tri.get_point(static_cast<size_t>(NorgIndex)));
354  }
355  }
356  }
357  else
358  {
359  if (_bc->GetBoundaryType() == Rectengular)
360  {
361  for (int i = static_cast<int>(counter); i < Tri.GetCorSize(); ++i)
362  {
363  int NorgIndex = Tri.GetOrgIndex(i);
364  if (NorgIndex < Nextra)
365  {
366  Vector2D norm = normalize(Tri.get_point(static_cast<size_t>(i)) -
367  Tri.get_point(static_cast<size_t>(NorgIndex)));
368  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + Tri.get_point(static_cast<size_t>(i))
369  - Tri.get_point(static_cast<size_t>(NorgIndex)) - 2 * norm * ScalarProd(norm,
370  CM[static_cast<size_t>(NorgIndex)] - Tri.get_point(static_cast<size_t>(NorgIndex)));
371  }
372  }
373  }
374  else // Half periodic case
375  {
376  double dx = _bc->GetGridBoundary(Right) - _bc->GetGridBoundary(Left);
377  for (int i = static_cast<int>(counter); i < Tri.GetCorSize(); ++i)
378  {
379  int NorgIndex = Tri.GetOrgIndex(i);
380  if (NorgIndex < Nextra)
381  {
382  double dx_temp = fastabs(Tri.get_point(static_cast<size_t>(i)) - Tri.get_point(static_cast<size_t>(NorgIndex)));
383  if (dx_temp<1.0001*dx && dx_temp*1.0001>dx)
384  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + (Tri.get_point(static_cast<size_t>(i)) -
385  Tri.get_point(static_cast<size_t>(NorgIndex)));
386  else
387  {
388  Vector2D norm = normalize(Tri.get_point(static_cast<size_t>(i)) - Tri.get_point(static_cast<size_t>(NorgIndex)));
389  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + Tri.get_point(static_cast<size_t>(i))
390  - Tri.get_point(static_cast<size_t>(NorgIndex)) - 2 * norm * ScalarProd(norm,
391  CM[static_cast<size_t>(NorgIndex)] - Tri.get_point(static_cast<size_t>(NorgIndex)));
392  }
393  }
394  }
395  }
396  }
397 }
398 
399 bool VoronoiMesh::legal_edge(Edge *e) //checks if both ends of the edge are outside the grid and that the edge doesn't cross the grid
400 {
401  bool n0 = e->neighbors.first < Tri.get_length();
402  bool n1 = e->neighbors.second < Tri.get_length();
403  if ( n0 || n1)
404  {
405  if (obc->GetBoundaryType() == Rectengular)
406  {
407  if (!n1 || !n0)
408  {
409 #ifdef RICH_MPI
410  if (!n0)
411  {
412  Vector2D p = Tri.get_point(n0);
413  if (p.x > obc->GetGridBoundary(Right) || p.x<obc->GetGridBoundary(Left) ||
414  p.y>obc->GetGridBoundary(Up) || p.y < obc->GetGridBoundary(Down))
415  {
416  if (GetOriginalIndex(n0) != GetOriginalIndex(n1))
417  return false;
418  }
419  }
420  else
421  {
422  Vector2D p = Tri.get_point(n1);
423  if (p.x > obc->GetGridBoundary(Right) || p.x<obc->GetGridBoundary(Left) ||
424  p.y>obc->GetGridBoundary(Up) || p.y < obc->GetGridBoundary(Down))
425  {
426  if (GetOriginalIndex(n0) != GetOriginalIndex(n1))
427  return false;
428  }
429  }
430 #else
431  if (GetOriginalIndex(e->neighbors.first) != GetOriginalIndex(e->neighbors.second))
432  return false;
433 #endif
434  }
435  }
436  return true;
437  }
438  else
439  return false;
440 }
441 
442 double VoronoiMesh::GetWidth(int index)const
443 {
444  return sqrt(GetVolume(index) / M_PI);
445 }
446 
447 vector<int> const& VoronoiMesh::GetCellEdges(int index) const
448 {
449  return mesh_vertices.at(static_cast<size_t>(index));
450 }
451 
452 double VoronoiMesh::GetVolume(int index) const
453 {
454  const Vector2D center = Tri.get_point(static_cast<size_t>(index));
455  double area = 0;
456  const size_t nloop = mesh_vertices[static_cast<size_t>(index)].size();
457  for (size_t i = 0; i < nloop; ++i)
458  {
459  const Vector2D p1 = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][static_cast<size_t>(i)])].vertices.first -
460  center;
461  const Vector2D p2 = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][static_cast<size_t>(i)])].vertices.second -
462  center;
463  area += 0.5*abs(ScalarProd(p1, zcross(p2)));
464  }
465  return area;
466 }
467 
468 Vector2D VoronoiMesh::CalcCellCM(size_t index) const
469 {
470  const Vector2D center = edges[static_cast<size_t>(mesh_vertices[index].front())].vertices.first;
471  Vector2D pc(0, 0);
472  double area = 0;
473  for (size_t i = 1; i < mesh_vertices[index].size(); i++)
474  {
475  const Edge& edge = edges[static_cast<size_t>(mesh_vertices[index][i])];
476  const Vector2D p1 = edge.vertices.first - center;
477  const Vector2D p2 = edge.vertices.second - center;
478  const double area_temp = 0.5*abs(ScalarProd(p1, zcross(p2)));
479  area += area_temp;
480  pc += (area_temp / 3.)*(center + edge.vertices.first + edge.vertices.second);
481  }
482  return pc / area;
483 }
484 
485 vector<Vector2D>& VoronoiMesh::GetMeshPoints(void)
486 {
487  return Tri.GetMeshPoints();
488 }
489 
490 vector<int> VoronoiMesh::Update(const vector<Vector2D>& pv, bool reorder)
491 {
492  // Clean_up last step
493  edges.clear();
494  GhostPoints.clear();
495  GhostProcs.clear();
496  vector<Vector2D> procpoints;
497  procpoints.push_back(Vector2D(obc->GetGridBoundary(Left), obc->GetGridBoundary(Down)));
498  procpoints.push_back(Vector2D(obc->GetGridBoundary(Right), obc->GetGridBoundary(Down)));
499  procpoints.push_back(Vector2D(obc->GetGridBoundary(Right), obc->GetGridBoundary(Up)));
500  procpoints.push_back(Vector2D(obc->GetGridBoundary(Left), obc->GetGridBoundary(Up)));
501 
502  vector<Vector2D> points = UpdatePoints(pv, obc);
503  // Check that points are all isnide domain
504  std::vector<std::pair<Vector2D, Vector2D> > chull = calc_procpoints(*obc);
505  CheckInput(chull, points);
506 
507  vector<int> HilbertIndeces;
508  if (reorder)
509  {
510  HilbertIndeces = HilbertOrder(points, static_cast<int>(pv.size()));
511  points = VectorValues(points, HilbertIndeces);
512  }
513 
514  Tri.update(points, chull);
515 
516  Nextra = static_cast<int>(Tri.ChangeCor().size());
517  vector<Edge> box_edges = obc->GetBoxEdges();
518  Tri.BuildBoundary(obc, box_edges);
519 
520  eps = 1e-8;
521  edges.clear();
522  GhostPoints.clear();
523  GhostProcs.clear();
524  build_v();
525 
526  if (logger)
527  logger->output(*this);
528 
529  CM.resize(Tri.getCor().size());
530  for (size_t i = 0; i < points.size(); ++i)
531  CM[i] = CalcCellCM(i);
532 
533  size_t counter = pv.size() + 3;
534  if (obc->GetBoundaryType() == Periodic)
535  {
536  for (int i = static_cast<int>(counter); i < Tri.GetCorSize(); ++i)
537  {
538  int NorgIndex = Tri.GetOrgIndex(i);
539  if (NorgIndex < Nextra)
540  {
541  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + (Tri.get_point(static_cast<size_t>(i)) -
542  Tri.get_point(static_cast<size_t>(NorgIndex)));
543  }
544  }
545  }
546  else
547  {
548  if (obc->GetBoundaryType() == Rectengular)
549  {
550  for (int i = static_cast<int>(counter); i < Tri.GetCorSize(); ++i)
551  {
552  int NorgIndex = Tri.GetOrgIndex(i);
553  if (NorgIndex < Nextra)
554  {
555  Vector2D norm = normalize(Tri.get_point(static_cast<size_t>(i)) - Tri.get_point(static_cast<size_t>(NorgIndex)));
556  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + Tri.get_point(static_cast<size_t>(i))
557  - Tri.get_point(static_cast<size_t>(NorgIndex)) - 2 * norm * ScalarProd(norm,
558  CM[static_cast<size_t>(NorgIndex)] - Tri.get_point(static_cast<size_t>(NorgIndex)));
559  }
560  }
561  }
562  else // Half periodic case
563  {
564  double dx = obc->GetGridBoundary(Right) - obc->GetGridBoundary(Left);
565  for (int i = static_cast<int>(counter); i < Tri.GetCorSize(); ++i)
566  {
567  int NorgIndex = Tri.GetOrgIndex(i);
568  if (NorgIndex < Nextra)
569  {
570  double dx_temp = fastabs(Tri.get_point(static_cast<size_t>(i)) - Tri.get_point(static_cast<size_t>(NorgIndex)));
571  if (dx_temp<1.0001*dx && dx_temp*1.0001>dx)
572  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + (Tri.get_point(static_cast<size_t>(i)) -
573  Tri.get_point(static_cast<size_t>(NorgIndex)));
574  else
575  {
576  Vector2D norm = normalize(Tri.get_point(static_cast<size_t>(i)) - Tri.get_point(static_cast<size_t>(NorgIndex)));
577  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + Tri.get_point(static_cast<size_t>(i))
578  - Tri.get_point(static_cast<size_t>(NorgIndex)) - 2 * norm * ScalarProd(norm,
579  CM[static_cast<size_t>(NorgIndex)] - Tri.get_point(static_cast<size_t>(NorgIndex)));
580  }
581  }
582  }
583  }
584  }
585  return HilbertIndeces;
586 }
587 
588 vector<int> VoronoiMesh::GetNeighbors(int index)const
589 {
590  vector<int> res(mesh_vertices[static_cast<size_t>(index)].size());
591  for (size_t i = 0; i < res.size(); ++i)
592  res[i] = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][i])].neighbors.first != index ?
593  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][i])].neighbors.first :
594  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][i])].neighbors.second;
595  return res;
596 }
597 
598 void VoronoiMesh::GetNeighbors(int index, vector<int> &neigh)const
599 {
600  neigh.resize(mesh_vertices[static_cast<size_t>(index)].size());
601  size_t N = neigh.size();
602  for (size_t i = 0; i < N; ++i)
603  neigh[i] = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][i])].neighbors.first != index ?
604  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][i])].neighbors.first :
605  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][i])].neighbors.second;
606 }
607 
608 
609 vector<int> VoronoiMesh::GetLiteralNeighbors(int index)const
610 {
611  int n = static_cast<int>(mesh_vertices[static_cast<size_t>(index)].size());
612  vector<int> res;
613  res.reserve(static_cast<size_t>(n));
614  for (int i = 0; i < n; ++i)
615  {
616  int other = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][static_cast<size_t>(i)])].neighbors.first;
617  if (other != index)
618  {
619  if (other > -1)
620  res.push_back(other);
621  }
622  else
623  {
624  if (other > -1)
625  other = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(index)][static_cast<size_t>(i)])].neighbors.second;
626  res.push_back(other);
627  }
628  }
629  return res;
630 }
631 
632 VoronoiMesh* VoronoiMesh::clone(void)const
633 {
634  return new VoronoiMesh(*this);
635 }
636 
637 int VoronoiMesh::GetPointNo(void) const
638 {
639  return Tri.get_length();
640 }
641 
642 void VoronoiMesh::SetPointNo(int N)
643 {
644  Tri.Changelength(N);
645  return;
646 }
647 
648 Vector2D VoronoiMesh::GetMeshPoint(int index) const
649 {
650  return Tri.get_point(static_cast<size_t>(index));
651 }
652 
653 int VoronoiMesh::GetTotalSidesNumber(void) const
654 {
655  return static_cast<int>(edges.size());
656 }
657 
658 const vector<Edge>& VoronoiMesh::getAllEdges(void) const
659 {
660  return edges;
661 }
662 
663 Edge const& VoronoiMesh::GetEdge(int index) const
664 {
665  return edges[static_cast<size_t>(index)];
666 }
667 
668 Vector2D const& VoronoiMesh::GetCellCM(int index) const
669 {
670  return CM[static_cast<size_t>(index)];
671 }
672 
673 vector<Vector2D>& VoronoiMesh::GetAllCM(void)
674 {
675  return CM;
676 }
677 
678 void VoronoiMesh::FindIntersectingOuterPoints(vector<Edge> const&box_edges, vector<vector<int> >
679  &boxduplicate, vector<vector<int> > const&firstduplicated)
680 {
681  int n = static_cast<int>(box_edges.size());
682  boxduplicate.resize(static_cast<size_t>(n));
683  int N = static_cast<int>(mesh_vertices.size());
684  if (N < 20)
685  {
686  for (int i = 0; i < n; ++i)
687  for (int j = 0; j < N; ++j)
688  boxduplicate[static_cast<size_t>(i)].push_back(j);
689  return;
690  }
691 
692  N = static_cast<int>(firstduplicated.size());
693  for (int i = 0; i < N; ++i)
694  {
695  n = static_cast<int>(firstduplicated[static_cast<size_t>(i)].size());
696  for (int j = 0; j < n; ++j)
697  {
698  vector<int> temp = CellIntersectOuterBoundary(box_edges, firstduplicated[static_cast<size_t>(i)][static_cast<size_t>(j)]);
699  int jj = static_cast<int>(temp.size());
700  if (jj > 0)
701  {
702  for (int k = 0; k < jj; ++k)
703  boxduplicate[static_cast<size_t>(temp[static_cast<size_t>(k)])].push_back(firstduplicated[static_cast<size_t>(i)][static_cast<size_t>(j)]);
704  }
705  }
706  }
707  n = static_cast<int>(box_edges.size());
708  for (int i = 0; i < n; ++i)
709  {
710  if (!boxduplicate[static_cast<size_t>(i)].empty())
711  {
712  sort(boxduplicate[static_cast<size_t>(i)].begin(), boxduplicate[static_cast<size_t>(i)].end());
713  boxduplicate[static_cast<size_t>(i)] = unique(boxduplicate[static_cast<size_t>(i)]);
714  }
715  }
716 }
717 
718 void VoronoiMesh::FindIntersectingPoints(vector<Edge> const &box_edges,
719  vector<vector<int> > &toduplicate)
720 {
721  int n = static_cast<int>(box_edges.size());
722  toduplicate.resize(static_cast<size_t>(n));
723  int N = static_cast<int>(mesh_vertices.size());
724  if (N < 20)
725  {
726  for (int i = 0; i < n; ++i)
727  for (int j = 0; j < N; ++j)
728  toduplicate[static_cast<size_t>(i)].push_back(j);
729  return;
730  }
731 
732  for (int i = 0; i < N; ++i)
733  {
734  vector<int> temp;
735  temp = CellIntersectBoundary(box_edges, i);
736  int j = static_cast<int>(temp.size());
737  if (j > 0)
738  {
739  for (int k = 0; k < j; ++k)
740  toduplicate[static_cast<size_t>(temp[static_cast<size_t>(k)])].push_back(i);
741  }
742  }
743  for (int i = 0; i < n; ++i)
744  {
745  if (!toduplicate[static_cast<size_t>(i)].empty())
746  {
747  sort(toduplicate[static_cast<size_t>(i)].begin(), toduplicate[static_cast<size_t>(i)].end());
748  toduplicate[static_cast<size_t>(i)] = unique(toduplicate[static_cast<size_t>(i)]);
749  }
750  }
751 }
752 
753 vector<int> VoronoiMesh::CellIntersectBoundary(vector<Edge> const&box_edges, int cell)
754 {
755  int ncell = static_cast<int>(mesh_vertices[static_cast<size_t>(cell)].size());
756  int nbox = static_cast<int>(box_edges.size());
757  vector<int> res;
758  Vector2D intersect;
759  for (int i = 0; i < ncell; ++i)
760  {
761  for (int j = 0; j < nbox; ++j)
762  {
763  if (SegmentIntersection(box_edges[static_cast<size_t>(j)], edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(cell)][static_cast<size_t>(i)])],
764  intersect))
765  res.push_back(j);
766  }
767  }
768  sort(res.begin(), res.end());
769  res = unique(res);
770  int nintersect = static_cast<int>(res.size());
771  if (nintersect > 1)
772  {
773  vector<Vector2D> cpoints;
774  ConvexHull(cpoints, *this, cell);
775  for (int i = 0; i < nbox; ++i)
776  if (PointInCell(cpoints, box_edges[static_cast<size_t>(i)].vertices.first) ||
777  PointInCell(cpoints, box_edges[static_cast<size_t>(i)].vertices.second))
778  res.push_back(i);
779  sort(res.begin(), res.end());
780  res = unique(res);
781  }
782  return res;
783 }
784 
785 vector<int> VoronoiMesh::CellIntersectOuterBoundary(vector<Edge> const&box_edges, int cell)
786 {
787  int ncell = static_cast<int>(mesh_vertices[static_cast<size_t>(cell)].size());
788  int nbox = static_cast<int>(box_edges.size());
789  vector<int> res;
790  Vector2D intersect;
791  boost::array<Vector2D, 3> tocheck;
792  for (int i = 0; i < ncell; ++i)
793  {
794  for (int j = 0; j < nbox; ++j)
795  {
796  if (SegmentIntersection(box_edges[static_cast<size_t>(j)], edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(cell)][static_cast<size_t>(i)])],
797  intersect))
798  {
799  double r = sqrt(edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(cell)][static_cast<size_t>(i)])].GetLength()*
800  box_edges[static_cast<size_t>(j)].GetLength());
801  double eps1 = 1e-7;
802  if (abs(orient2d(TripleConstRef<Vector2D>
803  (box_edges[static_cast<size_t>(j)].vertices.second - box_edges[static_cast<size_t>(j)].vertices.first,
804  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(cell)][static_cast<size_t>(i)])].vertices.second - box_edges[static_cast<size_t>(j)].vertices.first,
805  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(cell)][static_cast<size_t>(i)])].vertices.first - box_edges[static_cast<size_t>(j)].vertices.first)))
806  < r*r*eps1)
807  continue;
808  if (DistanceToEdge(edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(cell)][static_cast<size_t>(i)])].vertices.first,
809  box_edges[static_cast<size_t>(j)]) < eps1*r)
810  continue;
811  if (DistanceToEdge(edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(cell)][static_cast<size_t>(i)])].vertices.second,
812  box_edges[static_cast<size_t>(j)]) < eps1*r)
813  continue;
814  if ((intersect.distance(edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(cell)][static_cast<size_t>(i)])].vertices.first)
815  > eps1*r) && (intersect.distance(edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(cell)][static_cast<size_t>(i)])].vertices.second)
816  > eps1*r))
817  res.push_back(j);
818  }
819  }
820  }
821  sort(res.begin(), res.end());
822  res = unique(res);
823  return res;
824 }
825 
826 
827 bool VoronoiMesh::CloseToBorder(int point, int &border)
828 {
829  int olength = Tri.GetOriginalLength();
830  int n = static_cast<int>(mesh_vertices[static_cast<size_t>(point)].size());
831  for (int i = 0; i < n; ++i)
832  {
833  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].neighbors.second == point)
834  border = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].neighbors.first;
835  else
836  border = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].neighbors.second;
837  if (border > olength)
838  return true;
839  }
840  return false;
841 }
842 
843 vector<int> VoronoiMesh::GetBorderingCells(vector<int> const& copied,
844  vector<int> const& totest, int tocheck, vector<int> tempresult, int outer)
845 {
846  int border, test;
847  int olength = Tri.GetOriginalLength();
848  tempresult.push_back(tocheck);
849  sort(tempresult.begin(), tempresult.end());
850  int n = static_cast<int>(mesh_vertices[static_cast<size_t>(tocheck)].size());
851  for (int i = 0; i < n; ++i)
852  {
853  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(tocheck)][static_cast<size_t>(i)])].neighbors.second == tocheck)
854  test = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(tocheck)][static_cast<size_t>(i)])].neighbors.first;
855  else
856  test = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(tocheck)][static_cast<size_t>(i)])].neighbors.second;
857  if (test >= olength)
858  continue;
859  if (test < 0)
860  continue;
861  if (CloseToBorder(test, border))
862  if (border == outer)
863  if (!binary_search(copied.begin(), copied.end(), test) &&
864  !binary_search(totest.begin(), totest.end(), test) &&
865  !binary_search(tempresult.begin(), tempresult.end(), test))
866  tempresult = GetBorderingCells(copied, totest, test, tempresult, outer);
867  }
868  return tempresult;
869 }
870 
871 void VoronoiMesh::GetAdditionalBoundary(vector<vector<int> > &copied,
872  vector<vector<int> > &neighbors, vector<vector<int> > &totest)
873 {
874  int nsides = static_cast<int>(copied.size());
875  // Get all the neighbors
876  neighbors.clear();
877  neighbors.resize(static_cast<size_t>(nsides));
878  for (int i = 0; i < nsides; ++i)
879  {
880  sort(copied[static_cast<size_t>(i)].begin(), copied[static_cast<size_t>(i)].end());
881  // look if there are boundary points neighbors
882  int n = static_cast<int>(totest[static_cast<size_t>(i)].size());
883  for (int j = 0; j < n; ++j)
884  {
885  if (totest[static_cast<size_t>(i)][static_cast<size_t>(j)] == -1)
886  continue;
887  vector<int> toadd;
888  int outer = 0;
889  if (CloseToBorder(totest[static_cast<size_t>(i)][static_cast<size_t>(j)], outer))
890  toadd = GetBorderingCells(copied[static_cast<size_t>(i)], totest[static_cast<size_t>(i)], totest[static_cast<size_t>(i)][static_cast<size_t>(j)], toadd, outer);
891  int nn = static_cast<int>(toadd.size());
892  for (int k = 0; k < nn; ++k)
893  neighbors[static_cast<size_t>(i)].push_back(toadd[static_cast<size_t>(k)]);
894  }
895  sort(neighbors[static_cast<size_t>(i)].begin(), neighbors[static_cast<size_t>(i)].end());
896  neighbors[static_cast<size_t>(i)] = unique(neighbors[static_cast<size_t>(i)]);
897  neighbors[static_cast<size_t>(i)] = RemoveList(neighbors[static_cast<size_t>(i)], copied[static_cast<size_t>(i)]);
898  }
899 }
900 
901 void VoronoiMesh::GetRealNeighbor(vector<int> &result, int point) const
902 {
903  result.reserve(7);
904  int n = static_cast<int>(mesh_vertices[static_cast<size_t>(point)].size());
905  int olength = Tri.GetOriginalLength();
906  for (int i = 0; i < n; ++i)
907  {
908  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].neighbors.first == point)
909  {
910  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].neighbors.second > -1 &&
911  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].neighbors.second < olength)
912  result.push_back(edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].neighbors.second);
913  }
914  else
915  {
916  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].neighbors.first > -1 &&
917  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].neighbors.first < olength)
918  result.push_back(edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].neighbors.first);
919  }
920  }
921  sort(result.begin(), result.end());
922  result = unique(result);
923 }
924 
925 void VoronoiMesh::GetNeighborNeighbors(vector<int> &result, int point) const
926 {
927  vector<int> neigh = GetNeighbors(point);
928  result.clear();
929  result.reserve(25);
930  //GetRealNeighbor(neigh,point);
931  int n = static_cast<int>(neigh.size());
932  for (int i = 0; i < n; ++i)
933  {
934 #ifdef RICH_MPI
935  if (GetOriginalIndex(neigh[static_cast<size_t>(i)]) > GetPointNo())
936  continue;
937 #endif
938  vector<int> temp = GetNeighbors(GetOriginalIndex(neigh[static_cast<size_t>(i)]));
939  for (size_t j = 0; j < temp.size(); ++j)
940  result.push_back(GetOriginalIndex(temp[j]));
941  }
942  sort(result.begin(), result.end());
943  result = unique(result);
944  // Remove self point
945  RemoveVal(result, point);
946 }
947 
948 void VoronoiMesh::GetNeighborNeighborsMPI(vector<int> &result, int point)
949 {
950  vector<int> neigh;
951  result.clear();
952  GetRealNeighbor(neigh, point);
953  int n = static_cast<int>(neigh.size());
954  for (int i = 0; i < n; ++i)
955  {
956  vector<int> temp;
957  GetRealNeighbor(temp, neigh[static_cast<size_t>(i)]);
958  int N = static_cast<int>(temp.size());
959  for (int j = 0; j < N; ++j)
960  result.push_back(temp[static_cast<size_t>(j)]);
961  }
962  sort(result.begin(), result.end());
963  unique(result);
964 }
965 
966 void VoronoiMesh::GetCorners(vector<vector<int> > &copied,
967  vector<vector<int> > &result)
968 {
969  // copied should be sorted already
970  int nsides = static_cast<int>(copied.size());
971  result.clear();
972  OrgCorner.clear();
973  OrgCorner.resize(static_cast<size_t>(nsides));
974  result.resize(static_cast<size_t>(nsides));
975  vector<vector<int> > toadd(static_cast<size_t>(nsides));
976  for (int i = 0; i < nsides; ++i)
977  {
978  int n = static_cast<int>(copied[static_cast<size_t>(i)].size());
979  for (int j = 0; j < n; ++j)
980  {
981  if (binary_search(copied[static_cast<size_t>((i + 1) % nsides)].begin(), copied[static_cast<size_t>((i + 1) % nsides)].end(),
982  copied[static_cast<size_t>(i)][static_cast<size_t>(j)]))
983  {
984  vector<int> temp;
985  GetNeighborNeighborsMPI(temp, copied[static_cast<size_t>(i)][static_cast<size_t>(j)]);
986  result[static_cast<size_t>(i)].insert(result[static_cast<size_t>(i)].end(), temp.begin(), temp.end());
987  temp = AddPointsAlongEdge(static_cast<size_t>(copied[static_cast<size_t>(i)][static_cast<size_t>(j)]), copied, i);
988  toadd[static_cast<size_t>((i + 1) % nsides)].insert(toadd[static_cast<size_t>((i + 1) % nsides)].end(), temp.begin(),
989  temp.end());
990  temp = GetNeighbors(copied[static_cast<size_t>(i)][static_cast<size_t>(j)]);
991  for (vector<int>::iterator it = temp.begin(); it != temp.end(); ++it)
992  if (*it<GetPointNo() && *it>-1)
993  OrgCorner[static_cast<size_t>(i)].push_back(*it);
994  OrgCorner[static_cast<size_t>(i)].push_back(copied[static_cast<size_t>(i)][static_cast<size_t>(j)]);
995  }
996  if (binary_search(copied[static_cast<size_t>((i - 1 + nsides) % nsides)].begin(), copied[static_cast<size_t>((i - 1 + nsides) % nsides)].end(),
997  copied[static_cast<size_t>(i)][static_cast<size_t>(j)]))
998  {
999  vector<int> temp;
1000  GetNeighborNeighborsMPI(temp, copied[static_cast<size_t>(i)][static_cast<size_t>(j)]);
1001  result[static_cast<size_t>((i - 1 + nsides) % nsides)].insert(result[static_cast<size_t>((i - 1 + nsides) % nsides)].end()
1002  , temp.begin(), temp.end());
1003  temp = AddPointsAlongEdge(static_cast<size_t>(copied[static_cast<size_t>(i)][static_cast<size_t>(j)]), copied, i);
1004  toadd[static_cast<size_t>((i - 1 + nsides) % nsides)].insert(toadd[static_cast<size_t>((i - 1 + nsides) % nsides)].end(),
1005  temp.begin(), temp.end());
1006  temp = GetNeighbors(copied[static_cast<size_t>(i)][static_cast<size_t>(j)]);
1007  for (vector<int>::iterator it = temp.begin(); it != temp.end(); ++it)
1008  if (*it<GetPointNo() && *it>-1)
1009  OrgCorner[static_cast<size_t>((i - 1 + nsides) % nsides)].push_back(*it);
1010  OrgCorner[static_cast<size_t>((i - 1 + nsides) % nsides)].push_back(copied[static_cast<size_t>(i)][static_cast<size_t>(j)]);
1011  }
1012  }
1013  }
1014  for (int i = 0; i < nsides; ++i)
1015  {
1016  copied[static_cast<size_t>(i)].insert(copied[static_cast<size_t>(i)].end(), toadd[static_cast<size_t>(i)].begin(), toadd[static_cast<size_t>(i)].end());
1017  sort(copied[static_cast<size_t>(i)].begin(), copied[static_cast<size_t>(i)].end());
1018  copied[static_cast<size_t>(i)] = unique(copied[static_cast<size_t>(i)]);
1019  sort(result[static_cast<size_t>(i)].begin(), result[static_cast<size_t>(i)].end());
1020  result[static_cast<size_t>(i)] = unique(result[static_cast<size_t>(i)]);
1021  if (!OrgCorner[static_cast<size_t>(i)].empty())
1022  {
1023  sort(OrgCorner[static_cast<size_t>(i)].begin(), OrgCorner[static_cast<size_t>(i)].end());
1024  OrgCorner[static_cast<size_t>(i)] = unique(OrgCorner[static_cast<size_t>(i)]);
1025  }
1026  }
1027 }
1028 
1029 void VoronoiMesh::GetToTest(vector<vector<int> > &copied, vector<vector<int> > &totest)
1030 {
1031  int nsides = static_cast<int>(copied.size());
1032  int olength = Tri.GetOriginalLength();
1033  // sort the vectors
1034  for (int i = 0; i < nsides; ++i)
1035  sort(copied[static_cast<size_t>(i)].begin(), copied[static_cast<size_t>(i)].end());
1036  totest.resize(static_cast<size_t>(nsides));
1037  int test = 0;
1038  for (int i = 0; i < nsides; ++i)
1039  {
1040  vector<int> totest2;
1041  int ncopy = static_cast<int>(copied[static_cast<size_t>(i)].size());
1042  for (int j = 0; j < ncopy; ++j)
1043  {
1044  int n = static_cast<int>(mesh_vertices[static_cast<size_t>(copied[static_cast<size_t>(i)][static_cast<size_t>(j)])].size());
1045  for (int k = 0; k < n; ++k)
1046  {
1047  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(copied[static_cast<size_t>(i)][static_cast<size_t>(j)])][static_cast<size_t>(k)])].neighbors.first ==
1048  copied[static_cast<size_t>(i)][static_cast<size_t>(j)])
1049  test = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(copied[static_cast<size_t>(i)][static_cast<size_t>(j)])][static_cast<size_t>(k)])].neighbors.second;
1050  else
1051  test = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(copied[static_cast<size_t>(i)][static_cast<size_t>(j)])][static_cast<size_t>(k)])].neighbors.first;
1052  if (test < olength)
1053  totest2.push_back(test);
1054  }
1055  }
1056  sort(totest2.begin(), totest2.end());
1057  totest2 = unique(totest2);
1058  totest[static_cast<size_t>(i)] = totest2;
1059  }
1060 }
1061 
1062 vector<int> VoronoiMesh::FindEdgeStartConvex(int point)
1063 {
1064  int n = static_cast<int>(mesh_vertices[static_cast<size_t>(point)].size());
1065  Vector2D min_point;
1066  int min_index = 0, p_index;
1067  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(0)])].vertices.first.x <
1068  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][0])].vertices.second.x)
1069  {
1070  min_point = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][0])].vertices.first;
1071  p_index = 0;
1072  }
1073  else
1074  {
1075  min_point = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][0])].vertices.second;
1076  p_index = 1;
1077  }
1078  for (int i = 1; i < n; ++i)
1079  {
1080  double R = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].GetLength();
1081  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].vertices.first.x < (min_point.x - R*eps))
1082  {
1083  min_point = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].vertices.first;
1084  min_index = i;
1085  p_index = 0;
1086  }
1087  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].vertices.second.x < (min_point.x - R*eps))
1088  {
1089  min_point = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].vertices.second;
1090  min_index = i;
1091  p_index = 1;
1092  }
1093  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].vertices.first.x < (min_point.x + R*eps) &&
1094  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].vertices.first.y < min_point.y)
1095  {
1096  min_point = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].vertices.first;
1097  min_index = i;
1098  p_index = 0;
1099  }
1100 
1101  if (edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].vertices.second.x < (min_point.x + R*eps) &&
1102  edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].vertices.second.y < min_point.y)
1103  {
1104  min_point = edges[static_cast<size_t>(mesh_vertices[static_cast<size_t>(point)][static_cast<size_t>(i)])].vertices.second;
1105  min_index = i;
1106  p_index = 1;
1107  }
1108  }
1109  vector<int> res(2);
1110  res[0] = min_index;
1111  res[1] = p_index;
1112  return res;
1113 }
1114 
1115 void VoronoiMesh::ConvexEdgeOrder(void)
1116 {
1117  int n = static_cast<int>(mesh_vertices.size());
1118  for (int i = 0; i < n; ++i)
1119  {
1120  double R = GetWidth(i);
1121  vector<int> min_index = FindEdgeStartConvex(i);
1122  int p_loc = min_index[1];
1123  int edge_loc = mesh_vertices[static_cast<size_t>(i)][static_cast<size_t>(min_index[0])];
1124  int nedges = static_cast<int>(mesh_vertices[static_cast<size_t>(i)].size());
1125  std::list<int> elist;
1126  for (int j = 0; j < nedges; ++j)
1127  {
1128  if (j != min_index[0])
1129  elist.push_back(mesh_vertices[static_cast<size_t>(i)][static_cast<size_t>(j)]);
1130  }
1131  vector<int> new_order;
1132  new_order.reserve(static_cast<size_t>(nedges));
1133  new_order.push_back(edge_loc);
1134  for (int j = 0; j < nedges; ++j)
1135  {
1136  if (j == min_index[0])
1137  continue;
1138  int nlist = static_cast<int>(elist.size());
1139  std::list<int>::iterator it = elist.begin();
1140  for (int k = 0; k < nlist; ++k)
1141  {
1142  double temp0 = pair_member(edges[static_cast<size_t>(edge_loc)].vertices, (p_loc + 1) % 2).distance(edges[static_cast<size_t>(*it)].vertices.first);
1143  if (temp0 < eps*R)
1144  {
1145  p_loc = 0;
1146  edge_loc = *it;
1147  elist.erase(it);
1148  new_order.push_back(edge_loc);
1149  break;
1150  }
1151  double temp1 = pair_member(edges[static_cast<size_t>(edge_loc)].vertices, (p_loc + 1) % 2).distance(
1152  edges[static_cast<size_t>(*it)].vertices.second);
1153  if (temp1 < eps*R)
1154  {
1155  p_loc = 1;
1156  edge_loc = *it;
1157  new_order.push_back(edge_loc);
1158  elist.erase(it);
1159  break;
1160  }
1161  ++it;
1162  }
1163  }
1164  mesh_vertices[static_cast<size_t>(i)] = new_order;
1165  }
1166 }
1167 
1168 vector<Edge>& VoronoiMesh::GetAllEdges(void)
1169 {
1170  return edges;
1171 }
1172 
1173 vector<vector<int> >& VoronoiMesh::GetDuplicatedPoints(void)
1174 {
1175  return GhostPoints;
1176 }
1177 
1178 vector<vector<int> >const& VoronoiMesh::GetDuplicatedPoints(void)const
1179 {
1180  return GhostPoints;
1181 }
1182 
1183 vector<vector<int> >& VoronoiMesh::GetGhostIndeces(void)
1184 {
1185  return NGhostReceived;
1186 }
1187 
1188 vector<vector<int> >const& VoronoiMesh::GetGhostIndeces(void)const
1189 {
1190  return NGhostReceived;
1191 }
1192 
1193 int VoronoiMesh::GetTotalPointNumber(void)const
1194 {
1195  return static_cast<int>(Tri.getCor().size());
1196 }
1197 
1198 vector<int> VoronoiMesh::GetDuplicatedProcs(void)const
1199 {
1200  return GhostProcs;
1201 }
1202 
1203 vector<int> VoronoiMesh::GetSentProcs(void)const
1204 {
1205  return SentProcs;
1206 }
1207 
1208 vector<vector<int> >const& VoronoiMesh::GetSentPoints(void)const
1209 {
1210  return SentPoints;
1211 }
1212 
1213 // cpoints must be convex hull, checks if vec is inside cpoints
1214 bool PointInCell(vector<Vector2D> const& cpoints, Vector2D const& vec)
1215 {
1216  double Rmax = 0.0;
1217  double eps = 1e-8;
1218  size_t endp = cpoints.size();
1219  for (size_t i = 0; i < endp; ++i)
1220  Rmax = std::max(Rmax, fastabs(cpoints[i] - cpoints[(i + 1) % endp]));
1221  for (size_t i = 0; i < endp; ++i)
1222  {
1223  if (fastabs(cpoints[i] - cpoints[(i + 1) % endp]) < eps * Rmax)
1224  continue;
1225  if (orient2d(TripleConstRef<Vector2D>(cpoints[i],
1226  cpoints[(i + 1) % endp],
1227  vec)) < -0.0)
1228  return false;
1229  }
1230  return true;
1231 }
1232 
1233 // cpoints must be convex hull, checks if vec is inside cpoints
1234 bool PointInCell(std::vector<std::pair<Vector2D, Vector2D> > const& cpoints, Vector2D const& vec)
1235 {
1236  size_t endp = cpoints.size();
1237  for (size_t i = 0; i < endp; ++i)
1238  {
1239  if (orient2d(TripleConstRef<Vector2D>(cpoints[i].first,
1240  cpoints[i].second, vec)) < -0.0)
1241  return false;
1242  }
1243  return true;
1244 }
1245 
1246 // result is : minx, maxx, miny, maxy
1247 boost::array<double, 4> VoronoiMesh::FindMaxCellEdges(void)
1248 {
1249  int n = static_cast<int>(cell_edges.size());
1250  boost::array<double, 4> res;
1251  res[0] = min(cell_edges[0].vertices.first.x, cell_edges[0].vertices.second.x);
1252  res[1] = max(cell_edges[0].vertices.first.x, cell_edges[0].vertices.second.x);
1253  res[2] = min(cell_edges[0].vertices.first.y, cell_edges[0].vertices.second.y);
1254  res[3] = max(cell_edges[0].vertices.first.y, cell_edges[0].vertices.second.y);
1255  for (int i = 1; i < n; ++i)
1256  {
1257  res[0] = min(min(cell_edges[static_cast<size_t>(i)].vertices.first.x, cell_edges[static_cast<size_t>(i)].vertices.second.x), res[0]);
1258  res[1] = max(max(cell_edges[static_cast<size_t>(i)].vertices.first.x, cell_edges[static_cast<size_t>(i)].vertices.second.x), res[1]);
1259  res[2] = min(min(cell_edges[static_cast<size_t>(i)].vertices.first.y, cell_edges[static_cast<size_t>(i)].vertices.second.y), res[2]);
1260  res[3] = max(max(cell_edges[static_cast<size_t>(i)].vertices.first.y, cell_edges[static_cast<size_t>(i)].vertices.second.y), res[3]);
1261  }
1262  return res;
1263 }
1264 
1265 #ifdef RICH_MPI
1266 vector<int> VoronoiMesh::Update
1267 (const vector<Vector2D>& points,
1268  const Tessellation& vproc, bool reorder)
1269 {
1270  NGhostReceived.clear();
1271  int rank;
1272  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1273  vector<int> cedges;
1274  ConvexEdges(cedges, vproc, rank);
1275  cell_edges.clear();
1276  for (size_t i = 0; i < cedges.size(); ++i)
1277  cell_edges.push_back(vproc.GetEdge(cedges[i]));
1278 
1279  // Get the convex hull of the cell
1280  std::vector<std::pair<Vector2D, Vector2D> > cpoints;
1281  ConvexHull(cpoints, vproc, rank);
1282 
1283  vector<int> HilbertIndeces;
1284  vector<Vector2D> newcor;
1285  if (reorder)
1286  {
1287  HilbertIndeces = HilbertOrder(points, static_cast<int>(points.size()));
1288  newcor = VectorValues(points, HilbertIndeces);
1289  }
1290  else
1291  newcor = points;
1292 
1293 
1294  // Did points move between procs?
1295  newcor = UpdateMPIPoints(vproc, rank, newcor, obc, selfindex, SentProcs, SentPoints);
1296 
1297  //Build the delaunay
1298  try
1299  {
1300  Tri.build_delaunay(newcor, cpoints);
1301  }
1302  catch (UniversalError & eo)
1303  {
1304  size_t bad_index = static_cast<size_t>(eo.GetValues()[0]);
1305  eo.AddEntry("rank", rank);
1306  eo.AddEntry("original point in cor x", points[bad_index].x);
1307  eo.AddEntry("original point in cor y", points[bad_index].y);
1308  vector<int> edge_index = vproc.GetCellEdges(rank);
1309  for (size_t i = 0; i < edge_index.size(); ++i)
1310  {
1311  eo.AddEntry("edge number", static_cast<double>(edge_index[i]));
1312  eo.AddEntry("Edge v0x", vproc.GetEdge(edge_index[i]).vertices.first.x);
1313  eo.AddEntry("Edge v0y", vproc.GetEdge(edge_index[i]).vertices.first.y);
1314  eo.AddEntry("Edge v1x", vproc.GetEdge(edge_index[i]).vertices.second.x);
1315  eo.AddEntry("Edge v1y", vproc.GetEdge(edge_index[i]).vertices.second.y);
1316  }
1317  throw;
1318  }
1319  Nextra = static_cast<int>(Tri.ChangeCor().size());
1320  eps = 1e-8;
1321  edges.clear();
1322  GhostPoints.clear();
1323  GhostProcs.clear();
1324  NGhostReceived.clear();
1325  try
1326  {
1327  pair<vector<vector<int> >, vector<int> > ptemp = Tri.BuildBoundary(obc, vproc, NGhostReceived);
1328  GhostPoints = ptemp.first;
1329  GhostProcs = ptemp.second;
1330  }
1331  catch (UniversalError & eo)
1332  {
1333  std::vector<double> xtemp;
1334  for (size_t i = 0; i < newcor.size(); ++i)
1335  xtemp.push_back(newcor[i].x);
1336  write_vector(xtemp, "xcor_" + int2str(rank) + ".txt", 15);
1337  xtemp.clear();
1338  for (size_t i = 0; i < newcor.size(); ++i)
1339  xtemp.push_back(newcor[i].y);
1340  write_vector(xtemp, "ycor_" + int2str(rank) + ".txt", 15);
1341  throw;
1342  }
1343  build_v();
1344 
1345  if (logger)
1346  logger->output(*this);
1347 
1348  size_t n = static_cast<size_t>(GetPointNo());
1349  CM.resize(Tri.getCor().size());
1350  for (size_t i = 0; i < n; ++i)
1351  CM[i] = CalcCellCM(i);
1352 
1353  if (obc->GetBoundaryType() == Periodic)
1354  {
1355  for (int i = Nextra; i < Tri.GetCorSize(); ++i)
1356  {
1357  int NorgIndex = Tri.GetOrgIndex(i);
1358  if (NorgIndex < Nextra)
1359  {
1360  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + (Tri.get_point(i) -
1361  Tri.get_point(static_cast<size_t>(NorgIndex)));
1362  }
1363  }
1364  }
1365  else
1366  {
1367  if (obc->GetBoundaryType() == Rectengular)
1368  {
1369  for (int i = Nextra; i < Tri.GetCorSize(); ++i)
1370  {
1371  int NorgIndex = Tri.GetOrgIndex(i);
1372  if (NorgIndex < Nextra)
1373  {
1374  Vector2D norm = normalize(Tri.get_point(i) - Tri.get_point(static_cast<size_t>(NorgIndex)));
1375  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + Tri.get_point(i)
1376  - Tri.get_point(static_cast<size_t>(NorgIndex)) - 2 * norm * ScalarProd(norm,
1377  CM[static_cast<size_t>(NorgIndex)] - Tri.get_point(static_cast<size_t>(NorgIndex)));
1378  }
1379  }
1380  }
1381  else // Half periodic case
1382  throw(UniversalError("No HalfPeriodic in MPI"));
1383  }
1384  // communicate the ghost CM
1385  vector<vector<Vector2D> > incoming = MPI_exchange_data(GhostProcs, GhostPoints, CM);
1386  // Add the recieved CM
1387  for (size_t i = 0; i < incoming.size(); ++i)
1388  for (size_t j = 0; j < incoming.at(i).size(); ++j)
1389  CM[NGhostReceived.at(i).at(j)] = incoming[i][j];
1390  // Fix for periodic changes
1391  if (obc->GetBoundaryType() == Periodic)
1392  {
1393  std::vector<Vector2D> & localcor = Tri.ChangeCor();
1394  incoming = MPI_exchange_data(GhostProcs, GhostPoints, localcor);
1395  for (size_t i = 0; i < incoming.size(); ++i)
1396  for (size_t j = 0; j < incoming.at(i).size(); ++j)
1397  CM[NGhostReceived.at(i).at(j)] += localcor[NGhostReceived.at(i).at(j)] - incoming[i][j];
1398  }
1399 
1400  return HilbertIndeces;
1401 }
1402 
1403 void VoronoiMesh::Initialise
1404 (vector<Vector2D> const& pv,
1405  Tessellation const& vproc,
1406  OuterBoundary const* outer, bool reorder)
1407 {
1408  NGhostReceived.clear();
1409  int rank;
1410  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1411  obc = outer;
1412  vector<int> cedges;
1413  ConvexEdges(cedges, vproc, rank);
1414  cell_edges.clear();
1415  for (size_t i = 0; i < cedges.size(); ++i)
1416  cell_edges.push_back(vproc.GetEdge(cedges[i]));
1417 
1418  vector<Vector2D> points;
1419  if (reorder)
1420  points = VectorValues(pv, HilbertOrder(pv, static_cast<int>(pv.size())));
1421  else
1422  points = pv;
1423  // Get the convex hull of the cell
1424  std::vector<std::pair<Vector2D, Vector2D> > cpoints;
1425  ConvexHull(cpoints, vproc, rank);
1426  // Check input
1427  CheckInput(cpoints, points);
1428 
1429  //Build the delaunay
1430  Tri.build_delaunay(points, cpoints);
1431  Nextra = static_cast<int>(Tri.ChangeCor().size());
1432  eps = 1e-8;
1433  edges.clear();
1434  GhostPoints.clear();
1435  GhostProcs.clear();
1436  NGhostReceived.clear();
1437  selfindex.resize(points.size());
1438  size_t npoints = points.size();
1439  for (size_t i = 0; i < npoints; ++i)
1440  selfindex[i] = i;
1441  try
1442  {
1443  pair<vector<vector<int> >, vector<int> > ptemp = Tri.BuildBoundary(obc, vproc, NGhostReceived);
1444  GhostPoints = ptemp.first;
1445  GhostProcs = ptemp.second;
1446  }
1447  catch (UniversalError & eo)
1448  {
1449  std::vector<double> xtemp;
1450  for (size_t i = 0; i < points.size(); ++i)
1451  xtemp.push_back(points[i].x);
1452  write_vector(xtemp, "xcor_" + int2str(rank) + ".txt", 15);
1453  xtemp.clear();
1454  for (size_t i = 0; i < points.size(); ++i)
1455  xtemp.push_back(points[i].y);
1456  write_vector(xtemp, "ycor_" + int2str(rank) + ".txt", 15);
1457  throw;
1458  }
1459  /*
1460  if(get_rank()==0){
1461  {
1462  ofstream dump("GhostProcs.txt");
1463  for(size_t i=0;i<GhostProcs.size();++i)
1464  dump << GhostProcs.at(i) << endl;
1465  dump.close();
1466  }
1467  {
1468  ofstream dump("GhostPoints.txt");
1469  for(size_t i=0;i<GhostPoints.size();++i){
1470  for(size_t j=0;j<GhostPoints.at(i).size();++j)
1471  dump << GhostPoints.at(i).at(j) << " ";
1472  dump << endl;
1473  }
1474  dump.close();
1475  }
1476  assert(false);
1477  }
1478  */
1479 
1480  build_v();
1481 
1482  if (logger)
1483  logger->output(*this);
1484 
1485  size_t n = static_cast<size_t>(GetPointNo());
1486  CM.resize(Tri.getCor().size());
1487  for (size_t i = 0; i < n; ++i)
1488  CM[i] = CalcCellCM(i);
1489 
1490  if (obc->GetBoundaryType() == Periodic)
1491  {
1492  for (int i = Nextra; i < Tri.GetCorSize(); ++i)
1493  {
1494  int NorgIndex = Tri.GetOrgIndex(i);
1495  if (NorgIndex < Nextra)
1496  {
1497  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + (Tri.get_point(i) -
1498  Tri.get_point(static_cast<size_t>(NorgIndex)));
1499  }
1500  }
1501  }
1502  else
1503  {
1504  if (obc->GetBoundaryType() == Rectengular)
1505  {
1506  for (int i = Nextra; i < Tri.GetCorSize(); ++i)
1507  {
1508  int NorgIndex = Tri.GetOrgIndex(i);
1509  if (NorgIndex < Nextra)
1510  {
1511  Vector2D norm = normalize(Tri.get_point(i) - Tri.get_point(static_cast<size_t>(NorgIndex)));
1512  CM[static_cast<size_t>(i)] = CM[static_cast<size_t>(NorgIndex)] + Tri.get_point(i)
1513  - Tri.get_point(static_cast<size_t>(NorgIndex)) - 2 * norm * ScalarProd(norm,
1514  CM[static_cast<size_t>(NorgIndex)] - Tri.get_point(static_cast<size_t>(NorgIndex)));
1515  }
1516  }
1517  }
1518  else // Half periodic case
1519  throw(UniversalError("No HalfPeriodic in MPI"));
1520  }
1521  // communicate the ghost CM
1522  vector<vector<Vector2D> > incoming = MPI_exchange_data(GhostProcs, GhostPoints, CM);
1523  // Add the recieved CM
1524  for (size_t i = 0; i < incoming.size(); ++i)
1525  for (size_t j = 0; j < incoming.at(i).size(); ++j)
1526  CM[NGhostReceived.at(i).at(j)] = incoming[i][j];
1527 
1528  // Fix for periodic changes
1529  if (obc->GetBoundaryType() == Periodic)
1530  {
1531  std::vector<Vector2D> & localcor = Tri.ChangeCor();
1532  incoming = MPI_exchange_data(GhostProcs, GhostPoints, localcor);
1533  for (size_t i = 0; i < incoming.size(); ++i)
1534  for (size_t j = 0; j < incoming.at(i).size(); ++j)
1535  CM[NGhostReceived.at(i).at(j)] += localcor[NGhostReceived.at(i).at(j)] - incoming[i][j];
1536  }
1537 }
1538 
1539 
1540 vector<Vector2D> VoronoiMesh::UpdateMPIPoints(Tessellation const& vproc, int rank,
1541  vector<Vector2D> &points, OuterBoundary const* obc, vector<size_t> &selfindex,
1542  vector<int> &sentproc, vector<vector<int> > &sentpoints)
1543 {
1544  vector<Vector2D> res;
1545  res.reserve(points.size());
1546  selfindex.clear();
1547  size_t npoints = points.size();
1548  size_t nproc = static_cast<size_t>(vproc.GetPointNo());
1549  const double dx = obc->GetGridBoundary(Right) - obc->GetGridBoundary(Left);
1550  const double dy = obc->GetGridBoundary(Up) - obc->GetGridBoundary(Down);
1551  vector<std::pair<Vector2D, Vector2D> > cproc;
1552  ConvexHull(cproc, vproc, rank);
1553  vector<int> neighbors = vproc.GetNeighbors(rank);
1554  vector<int> realneigh;
1555  std::vector<size_t> neigh_keep;
1556  vector<vector<std::pair<Vector2D, Vector2D> > > neigh_chull;
1557  sentpoints.clear();
1558  sentproc.clear();
1559  bool periodic = obc->GetBoundaryType() == Periodic;
1560  for (size_t i = 0; i < neighbors.size(); ++i)
1561  {
1562  if (static_cast<size_t>(neighbors[i]) < nproc || periodic)
1563  {
1564  if (static_cast<size_t>(neighbors[i]) >= nproc)
1565  {
1566  int temp_neigh = vproc.GetOriginalIndex(neighbors[i]);
1567  if (temp_neigh == rank)
1568  continue;
1569  realneigh.push_back(temp_neigh);
1570  vector< std::pair<Vector2D, Vector2D> > temp;
1571  ConvexHull(temp, vproc, temp_neigh);
1572  Vector2D to_add_neigh = vproc.GetMeshPoint(neighbors[i]) - vproc.GetMeshPoint(temp_neigh);
1573  for (size_t j = 0; j < temp.size(); ++j)
1574  {
1575  temp[j].first += to_add_neigh;
1576  temp[j].second += to_add_neigh;
1577  }
1578  neigh_chull.push_back(temp);
1579  sentproc.push_back(temp_neigh);
1580  neigh_keep.push_back(i);
1581  }
1582  else
1583  {
1584  realneigh.push_back(neighbors[i]);
1585  vector< std::pair<Vector2D, Vector2D> > temp;
1586  ConvexHull(temp, vproc, neighbors[i]);
1587  neigh_chull.push_back(temp);
1588  sentproc.push_back(neighbors[i]);
1589  neigh_keep.push_back(i);
1590  }
1591  }
1592  }
1593  neighbors = VectorValues(neighbors, neigh_keep);
1594  std::sort(sentproc.begin(), sentproc.end());
1595  sentproc = unique(sentproc);
1596  sentpoints.resize(sentproc.size());
1597 
1598  for (size_t i = 0; i < npoints; ++i)
1599  {
1600  bool good = false;
1601  Vector2D temp = points[i];
1602  if (PointInCell(cproc, temp)) // Check own cpu
1603  {
1604  res.push_back(temp);
1605  selfindex.push_back(i);
1606  good = true;
1607  }
1608  if (!good)
1609  {
1610  for (size_t j = 0; j < realneigh.size(); ++j) // check cpu neighbors
1611  {
1612  if((!periodic && PointInCell(neigh_chull[j], temp)) || (periodic && PointInCell(neigh_chull[j], temp)))
1613  {
1614  good = true;
1615  if (periodic && static_cast<size_t>(neighbors[j]) >= nproc) // Do we need to move point?
1616  points[i] += vproc.GetMeshPoint(realneigh[j]) - vproc.GetMeshPoint(neighbors[j]);
1617  size_t index = std::find(sentproc.begin(), sentproc.end(), realneigh[j]) - sentproc.begin();
1618  assert(index < sentproc.size());
1619  sentpoints[index].push_back(static_cast<int>(i));
1620  break;
1621  }
1622  }
1623  }
1624  if (!good)
1625  {
1626  vector<std::pair<Vector2D, Vector2D> > cellpoints;
1627  for (size_t j = 0; j < nproc; ++j) // Search all cpus
1628  {
1629  if (std::find(realneigh.begin(), realneigh.end(), j) != realneigh.end() || j == static_cast<size_t>(rank))
1630  continue;
1631  ConvexHull(cellpoints, vproc, static_cast<int>(j));
1632  if ((!periodic && PointInCell(cellpoints, temp)) || (periodic && PointInCell(cellpoints, temp)))
1633  {
1634  good = true;
1635  size_t index = std::find(sentproc.begin(), sentproc.end(), j) - sentproc.begin();
1636  if (index >= sentproc.size())
1637  {
1638  sentproc.push_back(static_cast<int>(j));
1639  sentpoints.push_back(vector<int>(1, static_cast<int>(i)));
1640  }
1641  else
1642  sentpoints[index].push_back(static_cast<int>(i));
1643  break;
1644  }
1645  }
1646  }
1647  // If periodic check for all instances of duplication
1648  if (!good && obc->GetBoundaryType() == Periodic)
1649  {
1650  vector<std::pair<Vector2D, Vector2D> > cellpoints;
1651  for (size_t j = 0; j < nproc; ++j)
1652  {
1653  ConvexHull(cellpoints, vproc, static_cast<int>(j));
1654  // Create periodic instances of cpu
1655  std::vector<Vector2D> moved_point(8, temp);
1656  moved_point[0] += Vector2D(dx, 0);
1657  moved_point[1] += Vector2D(-dx, 0);
1658  moved_point[2] += Vector2D(dx, dy);
1659  moved_point[3] += Vector2D(dx, -dy);
1660  moved_point[4] += Vector2D(-dx, dy);
1661  moved_point[5] += Vector2D(-dx, -dy);
1662  moved_point[6] += Vector2D(0, dy);
1663  moved_point[7] += Vector2D(0, -dy);
1664  // Check if inside
1665  for (size_t k = 0; k < 8; ++k)
1666  {
1667  if (PointInCell(cellpoints, moved_point[k]))
1668  {
1669  good = true;
1670  points[i] = moved_point[k];
1671  if (static_cast<int>(j) == rank)
1672  {
1673  res.push_back(moved_point[k]);
1674  selfindex.push_back(i);
1675  break;
1676  }
1677  size_t index = std::find(sentproc.begin(), sentproc.end(), j) - sentproc.begin();
1678  if (index >= sentproc.size())
1679  {
1680  sentproc.push_back(static_cast<int>(j));
1681  sentpoints.push_back(vector<int>(1, static_cast<int>(i)));
1682  }
1683  else
1684  sentpoints[index].push_back(static_cast<int>(i));
1685  break;
1686  }
1687  }
1688  }
1689  }
1690  if (good)
1691  continue;
1692  HDF5Logger log("verror" + int2str(rank) + ".h5");
1693  log.output(vproc);
1694  UniversalError eo("Point is not inside any processor");
1695  eo.AddEntry("CPU rank", rank);
1696  eo.AddEntry("Point number", static_cast<double>(i));
1697  eo.AddEntry("Point x cor", points[i].x);
1698  eo.AddEntry("Point y cor", points[i].y);
1699  eo.AddEntry("cpu points", static_cast<double>(rank));
1700  for (size_t jj = 0; jj < cproc.size(); ++jj)
1701  {
1702  eo.AddEntry("x1", cproc[jj].first.x);
1703  eo.AddEntry("y1", cproc[jj].first.y);
1704  eo.AddEntry("x2", cproc[jj].second.x);
1705  eo.AddEntry("y2", cproc[jj].second.y);
1706  eo.AddEntry("orient result", orient2d(TripleConstRef<Vector2D>(cproc[jj].first,
1707  cproc[jj].second, points[i])));
1708  }
1709  for (size_t jj = 0; jj < neigh_chull.size(); ++jj)
1710  {
1711  eo.AddEntry("Cell number", static_cast<double>(realneigh[jj]));
1712  for (size_t kk = 0; kk < neigh_chull[jj].size(); ++kk)
1713  {
1714  eo.AddEntry("x1", neigh_chull[jj][kk].first.x);
1715  eo.AddEntry("y1", neigh_chull[jj][kk].first.y);
1716  eo.AddEntry("x2", neigh_chull[jj][kk].second.x);
1717  eo.AddEntry("y2", neigh_chull[jj][kk].second.y);
1718  eo.AddEntry("orient result", orient2d(TripleConstRef<Vector2D>(neigh_chull[jj][kk].first,
1719  neigh_chull[jj][kk].second, points[i])));
1720  }
1721  }
1722  throw eo;
1723  }
1724  // Send/Recv the points
1725  // Communication
1726  int wsize;
1727  MPI_Comm_size(MPI_COMM_WORLD, &wsize);
1728  vector<int> totalk(static_cast<size_t>(wsize), 0);
1729  vector<int> scounts(totalk.size(), 1);
1730  for (size_t i = 0; i < sentproc.size(); ++i)
1731  totalk[sentproc[i]] = 1;
1732  int nrecv;
1733  MPI_Reduce_scatter(&totalk[0], &nrecv, &scounts[0], MPI_INT, MPI_SUM,
1734  MPI_COMM_WORLD);
1735 
1736  vector<MPI_Request> req(sentproc.size());
1737  for (size_t i = 0; i < sentproc.size(); ++i)
1738  MPI_Isend(&wsize, 1, MPI_INT, sentproc[i], 3, MPI_COMM_WORLD, &req[i]);
1739  vector<int> talkwithme;
1740  for (int i = 0; i < nrecv; ++i)
1741  {
1742  MPI_Status status;
1743  MPI_Recv(&wsize, 1, MPI_INT, MPI_ANY_SOURCE, 3, MPI_COMM_WORLD, &status);
1744  talkwithme.push_back(status.MPI_SOURCE);
1745  }
1746  if (req.size() > 0)
1747  MPI_Waitall(static_cast<int>(req.size()), &req[0], MPI_STATUSES_IGNORE);
1748  MPI_Barrier(MPI_COMM_WORLD);
1749  for (size_t i = 0; i < talkwithme.size(); ++i)
1750  {
1751  if (std::find(sentproc.begin(), sentproc.end(), talkwithme[i]) == sentproc.end())
1752  {
1753  sentproc.push_back(talkwithme[i]);
1754  sentpoints.push_back(vector<int>());
1755  }
1756  }
1757  // Point exchange
1758  vector<vector<Vector2D> > incoming(sentproc.size());
1759  vector<vector<double> > tosend(sentproc.size());
1760  vector<double> torecv;
1761  double dtemp = 0;
1762  Vector2D vtemp;
1763  req.clear();
1764  req.resize(sentproc.size());
1765  vector<int> indeces;
1766  vector<Vector2D> cortemphilbert;
1767  for (size_t i = 0; i < sentproc.size(); ++i)
1768  {
1769  const int dest = sentproc.at(i);
1770  if (!sentpoints.at(i).empty())
1771  {
1772  cortemphilbert = VectorValues(points, sentpoints.at(i));
1773  indeces = HilbertOrder(cortemphilbert, static_cast<int>(cortemphilbert.size()));
1774  tosend[i] = list_serialize(VectorValues(cortemphilbert, indeces));
1775  sentpoints[i] = VectorValues(sentpoints[i], indeces);
1776  }
1777  if (tosend[i].empty())
1778  MPI_Isend(&dtemp, 1, MPI_DOUBLE, dest, 1, MPI_COMM_WORLD, &req[i]);
1779  else
1780  MPI_Isend(&tosend[i][0], static_cast<int>(tosend[i].size()), MPI_DOUBLE, dest, 0, MPI_COMM_WORLD, &req[i]);
1781  }
1782  for (size_t i = 0; i < sentproc.size(); ++i)
1783  {
1784  MPI_Status status;
1785  MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
1786  int count;
1787  MPI_Get_count(&status, MPI_DOUBLE, &count);
1788  torecv.resize(static_cast<size_t>(count));
1789  MPI_Recv(&torecv[0], count, MPI_DOUBLE, status.MPI_SOURCE, status.MPI_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
1790  if (status.MPI_TAG == 0)
1791  {
1792  size_t location = static_cast<size_t>(std::find(sentproc.begin(), sentproc.end(), status.MPI_SOURCE) -
1793  sentproc.begin());
1794  if (location >= sentproc.size())
1795  throw UniversalError("Bad location in mpi exchange");
1796  incoming[location] = list_unserialize(torecv, vtemp);
1797  }
1798  else
1799  {
1800  if (status.MPI_TAG != 1)
1801  throw UniversalError("Recv bad mpi tag");
1802  }
1803  }
1804  if (req.size() > 0)
1805  MPI_Waitall(static_cast<int>(req.size()), &req[0], MPI_STATUSES_IGNORE);
1806  MPI_Barrier(MPI_COMM_WORLD);
1807 
1808  // Combine the vectors
1809  for (size_t i = 0; i < incoming.size(); ++i)
1810  for (size_t j = 0; j < incoming[i].size(); ++j)
1811  res.push_back(incoming[i][j]);
1812  return res;
1813 }
1814 
1815 
1816 
1817 #endif
unique
vector< T > unique(vector< T > const &v)
Returns a vector containing only unique elements.
Definition: utils.hpp:376
VoronoiMesh
Voronoi tessellation class.
Definition: VoronoiMesh.hpp:29
Vector2D::distance
double distance(Vector2D const &v1) const
Caluclates the distance from the Vector to v1.
Definition: geometry.cpp:13
VoronoiMesh::GetSentProcs
vector< int > GetSentProcs(void) const
Returns the indeces of the processors with whom points where exchanged.
Definition: VoronoiMesh.cpp:1203
VoronoiMesh::GetOriginalIndex
int GetOriginalIndex(int point) const
Returns the original index of the duplicated point.
Definition: VoronoiMesh.cpp:173
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
VoronoiMesh::GetEdge
Edge const & GetEdge(int index) const
Returns edge (interface between cells)
Definition: VoronoiMesh.cpp:663
VoronoiMesh::GetDuplicatedPoints
vector< vector< int > > & GetDuplicatedPoints(void)
Returns the indeces of the points that where sent to other processors as ghost points.
Definition: VoronoiMesh.cpp:1173
sort_index
void sort_index(const vector< T > &arr, vector< int > &res)
Returns the indeces of a sort.
Definition: utils.hpp:500
SegmentIntersection
bool SegmentIntersection(Edge const &edge1, Edge const &edge2, Vector2D &Intersection, double eps=1e-8)
Calculates the intersection of two edges.
Definition: Edge.cpp:45
VoronoiMesh::GetTotalPointNumber
int GetTotalPointNumber(void) const
Returns the total number of points (including ghost)
Definition: VoronoiMesh.cpp:1193
VoronoiMesh::Update
vector< int > Update(const vector< Vector2D > &points, const Tessellation &vproc, bool reorder=false)
Updates the tessellation.
Definition: VoronoiMesh.cpp:1267
Tessellation
Abstract class for tessellation.
Definition: tessellation.hpp:21
VoronoiMesh::GetSelfPoint
vector< size_t > GetSelfPoint(void) const
Returns the indeces of the points that remain with the processor after the ne processor mesh is built...
Definition: VoronoiMesh.cpp:182
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
write_vector
void write_vector(vector< double > const &v, string const &fname, int prec=6)
Writes a list of numbers to a file.
Definition: simple_io.cpp:15
VoronoiMesh::UpdateMPIPoints
vector< Vector2D > UpdateMPIPoints(Tessellation const &vproc, int rank, vector< Vector2D > &points, OuterBoundary const *obc, vector< size_t > &selfindex, vector< int > &sentproc, vector< vector< int > > &sentpoints)
Communicate position of mesh generating points between processes.
Definition: VoronoiMesh.cpp:1540
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
PointInCell
bool PointInCell(vector< Vector2D > const &cpoints, Vector2D const &vec)
Checks if a point is inside a Voronoi cell.
Definition: VoronoiMesh.cpp:1214
Delaunay::get_length
int get_length(void) const
Returns the number of points.
Definition: Delaunay.cpp:724
UniversalError::GetValues
std::vector< double > const & GetValues(void) const
Returns entry values.
Definition: universal_error.cpp:37
VoronoiMesh::GetSentPoints
vector< vector< int > > const & GetSentPoints(void) const
Returns the indeces of the points that where sent to other processors.
Definition: VoronoiMesh.cpp:1208
VoronoiMesh::GetNeighbors
vector< int > GetNeighbors(int index) const
Returns the indeces of the neighbors.
Definition: VoronoiMesh.cpp:588
HDF5Logger::output
void output(const VoronoiMesh &v)
Outputs information from the Voronoi tessellation.
Definition: hdf5_logger.cpp:35
Edge
Interface between two cells.
Definition: Edge.hpp:13
Tessellation::GetMeshPoint
virtual Vector2D GetMeshPoint(int index) const =0
Returns Position of mesh generating point.
VoronoiMesh::Initialise
void Initialise(vector< Vector2D > const &points, Tessellation const &vproc, OuterBoundary const *outer, bool reorder=true)
Initialise the tessellation.
Definition: VoronoiMesh.cpp:1404
voronoi_loggers::VoronoiLogger::output
virtual void output(const VoronoiMesh &v)
Outputs information from the Voronoi tessellation.
Definition: voronoi_logger.cpp:10
VoronoiMesh::GetLiteralNeighbors
vector< int > GetLiteralNeighbors(int index) const
Returns the list of neighbors including ghost points.
Definition: VoronoiMesh.cpp:609
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
VoronoiMesh::GetDuplicatedProcs
vector< int > GetDuplicatedProcs(void) const
Returns the indeces of the processors with whom ghost points where exchanged.
Definition: VoronoiMesh.cpp:1198
int2str
std::string int2str(int n)
Converts an integer to a string.
Definition: int2str.cpp:6
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)
pair_member
T pair_member(const std::pair< T, T > &p, int index)
Selects a member of std::pair.
Definition: utils.hpp:660
Delaunay::BuildBoundary
void BuildBoundary(OuterBoundary const *obc, vector< Edge > const &edges)
Builds the boundary points.
Definition: Delaunay.cpp:1167
VoronoiMesh::GetWidth
double GetWidth(int index) const
Returns the effective width of a cell.
Definition: VoronoiMesh.cpp:442
ConvexHull
void ConvexHull(vector< Vector2D > &result, Tessellation const &tess, int index)
Returns the ConvexHull for a set of points.
Definition: ConvexHull.cpp:31
Delaunay::get_num_facet
int get_num_facet(void) const
Returns the number of facets.
Definition: Delaunay.cpp:719
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
VoronoiMesh::NearBoundary
bool NearBoundary(int index) const
Returns if the cell is adjacent to a boundary.
Definition: VoronoiMesh.cpp:157
VoronoiMesh::GetNeighborNeighbors
void GetNeighborNeighbors(vector< int > &result, int point) const
Retrieves vicarious neighbors.
Definition: VoronoiMesh.cpp:925
VoronoiMesh::GetAllEdges
vector< Edge > & GetAllEdges(void)
Returns a reference to a list of all edges.
Definition: VoronoiMesh.cpp:1168
VoronoiMesh::GetPointNo
int GetPointNo(void) const
Get Total number of mesh generating points.
Definition: VoronoiMesh.cpp:637
VoronoiMesh::GetCellCM
Vector2D const & GetCellCM(int index) const
Returns Position of Cell&#39;s CM.
Definition: VoronoiMesh.cpp:668
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
Edge::vertices
std::pair< Vector2D, Vector2D > vertices
Points at the ends of the edge.
Definition: Edge.hpp:18
VoronoiMesh::GetCellEdges
vector< int > const & GetCellEdges(int index) const
Returns a reference to a vector<int> containing the indexes of the edges related to a cell...
Definition: VoronoiMesh.cpp:447
VoronoiMesh::getAllEdges
const vector< Edge > & getAllEdges(void) const
Returns reference to the list of all edges.
Definition: VoronoiMesh.cpp:658
Delaunay::GetOrgIndex
int GetOrgIndex(int index) const
Retrieves the original index of a point (in case a point was duplicated)
Definition: Delaunay.cpp:1963
hdf5_logger.hpp
Method for dumping tessellation data to hdf5 file.
reorder
void reorder(vector< T > &v, vector< size_t > const &order)
Reorder a vector according to an index vector (obtained from the &#39;ordered&#39; function) ...
Definition: HilbertOrder3D_Utils.hpp:39
Delaunay::get_point
Vector2D get_point(size_t index) const
Returns a point.
Definition: Delaunay.cpp:704
min
double min(vector< double > const &v)
Returns the minimal term in a vector.
Definition: utils.cpp:44
normalize
Vector2D normalize(const Vector2D &v)
Normalized a vector.
Definition: geometry.cpp:158
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
VoronoiMesh::GetTotalSidesNumber
int GetTotalSidesNumber(void) const
Returns the total number of faces.
Definition: VoronoiMesh.cpp:653
VoronoiMesh::clone
VoronoiMesh * clone(void) const
Cloning function.
Definition: VoronoiMesh.cpp:632
Edge::neighbors
std::pair< int, int > neighbors
Neighboring cells.
Definition: Edge.hpp:21
DistanceToEdge
double DistanceToEdge(Vector2D const &point, Edge const &edge)
Calculates the distance of a point to an edge.
Definition: Edge.cpp:31
VoronoiMesh::SetPointNo
void SetPointNo(int N)
Set the number of points.
Definition: VoronoiMesh.cpp:642
VoronoiMesh::VoronoiMesh
VoronoiMesh(void)
Class default constructor.
Definition: VoronoiMesh.cpp:187
Delaunay::GetOriginalLength
int GetOriginalLength(void) const
Returns the original length of the points (without duplicated points)
Definition: Delaunay.cpp:739
VoronoiMesh.hpp
Voronoi tessellation with MPI option.
VoronoiMesh::GetGhostIndeces
vector< vector< int > > & GetGhostIndeces(void)
Returns the indeces of each ghost point in the vector of points that the tessellation holds...
Definition: VoronoiMesh.cpp:1183
abs
double abs(Vector3D const &v)
Norm of a vector.
Definition: Vector3D.cpp:44
VoronoiMesh::GetVolume
double GetVolume(int index) const
Returns the volume of a cell.
Definition: VoronoiMesh.cpp:452
VoronoiMesh::GetMeshPoints
vector< Vector2D > & GetMeshPoints(void)
Returns a reference to the point vector.
Definition: VoronoiMesh.cpp:485
UniversalError::AddEntry
void AddEntry(std::string const &field, double value)
Adds an entry to the list.
Definition: universal_error.cpp:20
HDF5Logger
Writes tessellation data to hdf5 format.
Definition: hdf5_logger.hpp:14
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
VoronoiMesh::logger
voronoi_loggers::VoronoiLogger * logger
Diagnostics method.
Definition: VoronoiMesh.hpp:169
Vector2D
2D Mathematical vector
Definition: geometry.hpp:15
RemoveList
vector< T > RemoveList(vector< T > const &v, vector< T > const &list)
Returns only elements from vector v which are not in vector list, assumes list is sorted...
Definition: utils.hpp:416
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
VoronoiMesh::GetMeshPoint
Vector2D GetMeshPoint(int index) const
Returns Position of mesh generating point.
Definition: VoronoiMesh.cpp:648
Edge::GetLength
double GetLength(void) const
Returns the length of the edge.
Definition: Edge.cpp:26
Tessellation::GetNeighbors
virtual vector< int > GetNeighbors(int index) const =0
Returns the indeces of the neighbors.
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
ConvexEdges
void ConvexEdges(vector< int > &result, Tessellation const &tess, int index)
Returns the ConvexHull of the edges of a cell.
Definition: ConvexHull.cpp:98
fastabs
double fastabs(Vector2D const &v)
Norm of a vector, less accurate.
Definition: geometry.cpp:24
VoronoiMesh::CalcFaceVelocity
Vector2D CalcFaceVelocity(Vector2D wl, Vector2D wr, Vector2D rL, Vector2D rR, Vector2D f) const
Calculates the velocity of a single edge.
Definition: VoronoiMesh.cpp:150
VoronoiMesh::GetAllCM
vector< Vector2D > & GetAllCM(void)
Returns the center of masses of the cells.
Definition: VoronoiMesh.cpp:673
Vector2D::x
double x
Component in the x direction.
Definition: geometry.hpp:45
zcross
Vector2D zcross(Vector2D const &v)
Cross product of a vector in x,y plane with a unit vector in the z direction.
Definition: geometry.cpp:145
Delaunay::ChangeCor
vector< Vector2D > & ChangeCor(void)
Allows to change the cor.
Definition: Delaunay.cpp:681
Generated on Thu Jun 4 2020 18:12:17 by   doxygen 1.8.13