VPTissue Reference Manual
EditableMesh.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2011-2016 Universiteit Antwerpen
3  *
4  * Licensed under the EUPL, Version 1.1 or as soon they will be approved by
5  * the European Commission - subsequent versions of the EUPL (the "Licence");
6  * You may not use this work except in compliance with the Licence.
7  * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
8  *
9  * Unless required by applicable law or agreed to in writing, software
10  * distributed under the Licence is distributed on an "AS IS" basis,
11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12  * See the Licence for the specific language governing
13  * permissions and limitations under the Licence.
14  */
20 #include "EditableMesh.h"
21 
22 #include "algo/CellDivider.h"
23 #include "bio/Cell.h"
24 #include "bio/CellAttributes.h"
25 #include "bio/Edge.h"
26 #include "bio/Mesh.h"
27 #include "bio/MeshCheck.h"
28 #include "bio/NeighborNodes.h"
29 #include "bio/Node.h"
30 #include "bio/PBMBuilder.h"
31 #include "bio/Wall.h"
32 #include "bio/WallAttributes.h"
33 #include "editor/PolygonUtils.h"
35 
36 #include <algorithm>
37 
38 namespace SimPT_Sim { class Mesh; }
39 
40 namespace SimPT_Editor {
41 
42 using namespace std;
43 using namespace SimPT_Sim;
44 using namespace SimPT_Sim::Container;
45 
46 EditableMesh::EditableMesh(const boost::property_tree::ptree& pt)
47  : m_mesh(PBMBuilder().Build(pt)), m_cells(m_mesh->GetCells()),
48  m_nodes(m_mesh->GetNodes()), m_walls(m_mesh->GetWalls()) {}
49 
51 
52 bool EditableMesh::CanSplitCell(Node* node1, Node* node2, Cell* cell) const
53 {
54  list<Cell*> owners;
55  for (auto nb : m_mesh->GetNodeOwningNeighbors(node1)) {
56  if (nb.GetCell()->GetIndex() != -1) {
57  owners.push_back(nb.GetCell());
58  }
59  }
60  for (auto nb : m_mesh->GetNodeOwningNeighbors(node2)) {
61  if (nb.GetCell()->GetIndex() != -1) {
62  owners.push_back(nb.GetCell());
63  }
64  }
65  owners.sort([](Cell* c1, Cell* c2) {return c1 < c2;});
66 
67  auto eql = [](Cell* c1, Cell* c2) {
68  return c1 == c2;
69  };
70 
71  list<Cell*> commonOwners;
72  auto owner = adjacent_find(owners.begin(), owners.end(), eql);
73  while (owner != owners.end()) {
74  commonOwners.push_back(*owner);
75  owner = adjacent_find(++owner, owners.end(), eql);
76  }
77 
78  if (commonOwners.empty()) {
79  return false;
80  }
81  else if (cell != nullptr) {
82  if (find(commonOwners.begin(), commonOwners.end(), cell) != commonOwners.end()) {
83  commonOwners.clear();
84  commonOwners.push_back(cell);
85  }
86  else {
87  return false;
88  }
89  }
90 
91  QLineF cut((*node1)[0], (*node1)[1], (*node2)[0], (*node2)[1]);
92  for (Cell* c : commonOwners) {
93  QPolygonF polygon;
94  for (Node* n : c->GetNodes()) {
95  polygon.push_back(QPointF((*n)[0], (*n)[1]));
96  }
97 
98  auto polygons = PolygonUtils::SlicePolygon(polygon, cut);
99  if (polygons.size() == 2) {
100  bool ok = true;
101  auto it = polygons.begin();
102  while (it != polygons.end() && ok) {
103  if (!it->contains(cut.p1()) || !it->contains(cut.p2())) {
104  ok = false;
105  }
106  for (const QPointF& p : *it) {
107  if (!polygon.contains(p)) {
108  ok = false;
109  }
110  }
111  it++;
112  }
113  if (ok) {
114  return true;
115  }
116  }
117  }
118 
119  return false;
120 }
121 
122 Cell* EditableMesh::CreateCell(Node* node1, Node* node2, const QPointF& newPoint)
123 {
124  assert(node1->IsAtBoundary() && node2->IsAtBoundary()
125  && "The two given nodes aren't located at the boundary of the mesh.");
126  assert(find(m_mesh->GetBoundaryPolygon()->GetNodes().begin(), m_mesh->GetBoundaryPolygon()->GetNodes().end(), node1) != m_mesh->GetBoundaryPolygon()->GetNodes().end()
127  && find(m_mesh->GetBoundaryPolygon()->GetNodes().begin(), m_mesh->GetBoundaryPolygon()->GetNodes().end(), node1) != m_mesh->GetBoundaryPolygon()->GetNodes().end()
128  && "The two given nodes aren't located at the boundary of the mesh.");
129 
130  QPolygonF boundaryPoints;
131  for (Node* n : m_mesh->GetBoundaryPolygon()->GetNodes()) {
132  boundaryPoints.append(QPointF((*n)[0], (*n)[1]));
133  }
134 
135  list<QLineF> excluded;
136  list<QLineF> included({QLineF(QPointF((*node1)[0], (*node1)[1]), newPoint), QLineF(QPointF((*node2)[0], (*node2)[1]), newPoint)});
137  if (boundaryPoints.contains(newPoint) || boundaryPoints.containsPoint(newPoint, Qt::OddEvenFill)
138  || !IsConsistent(excluded, included)) {
139 
140  return nullptr;
141  }
142  else {
143  // There are two possible polygons which can be the cell.
144  // The two ways you can go around the boundary polygon.
145  list<QPolygonF> polygons({QPolygonF(), QPolygonF()});
146  auto current = polygons.begin();
147 
148  QPointF point1((*node1)[0], (*node1)[1]);
149  QPointF point2((*node2)[0], (*node2)[1]);
150 
151  for (const QPointF& p : boundaryPoints) {
152  current->append(p);
153  if (p == point1 || p == point2) {
154  current->append(newPoint);
155  if (++current == polygons.end()) {
156  current = polygons.begin();
157  }
158  current->append(p);
159  }
160  }
161 
162  // The smallest one is the right one, since the largest one will form the boundary polygon.
163  QPolygonF polygon = polygons.front().intersected(polygons.back());
164  if (polygon.front() == polygon.back()) {
165  polygon.pop_back(); //So the polygon will become open.
166  }
167 
168  // Make the polygon counterclockwise (in the Qt coordinate system).
169  if (PolygonUtils::IsClockwise(polygon)) {
170  reverse(polygon.begin(), polygon.end());
171  }
172 
173  // After this, the nodes are in the order in which they will appear in the wall
174  // between the new cell and the boundary polygon.
175  rotate(polygon.begin(), find(polygon.begin(), polygon.end(), newPoint), polygon.end());
176  if (*next(polygon.begin()) == point1) {
177  swap(node1, node2);
178  swap(point1, point2);
179  }
180  else {
181  assert(*next(polygon.begin()) == point2
182  && "One of the two points should be next to the new point.");
183  }
184 
185  //Build the new node.
186  auto newNode = m_mesh->BuildNode({{newPoint.x(), newPoint.y(), 0}}, true);
187 
188  //Order old boundary nodes to have the first node (start of the wall) at the front.
189  auto oldBoundaryNodes = m_mesh->GetBoundaryPolygon()->GetNodes();
190  auto node_it = find(oldBoundaryNodes.begin(), oldBoundaryNodes.end(), node1);
191  rotate(oldBoundaryNodes.begin(), node_it, oldBoundaryNodes.end());
192 
193  //Create the new boundary nodes with the first node (start of the wall) at the front.
194  auto& newBoundaryNodes = m_mesh->GetBoundaryPolygon()->GetNodes();
195  node_it = find(newBoundaryNodes.begin(), newBoundaryNodes.end(), node2);
196  rotate(newBoundaryNodes.begin(), node_it, newBoundaryNodes.end());
197  node_it = find(newBoundaryNodes.begin(), newBoundaryNodes.end(), node1);
198  newBoundaryNodes.erase(next(node_it), newBoundaryNodes.end());
199  newBoundaryNodes.push_back(newNode); // may invalidate node_it, so compute it again
200  node_it = find(newBoundaryNodes.begin(), newBoundaryNodes.end(), node1);
201  rotate(newBoundaryNodes.begin(), node_it, newBoundaryNodes.end());
202  m_mesh->GetBoundaryPolygon()->SetGeoDirty();
203 
204  // Update the walls.
205  if (m_walls.size() == 0) {
206  assert(m_cells.size() == 1 && "There should be two cells in the mesh.");
207 
208  Cell* cell = m_cells.front();
209 
210  Wall* newWall = m_mesh->BuildWall(node2, node1, cell, m_mesh->GetBoundaryPolygon());
211  newWall->GetC1()->AddWall(newWall);
212  newWall->GetC2()->AddWall(newWall);
213  m_mesh->UpdateNodeOwningWalls(newWall);
214 
215  newWall = m_mesh->BuildWall(node1, node2, cell, m_mesh->GetBoundaryPolygon());
216  newWall->GetC1()->AddWall(newWall);
217  newWall->GetC2()->AddWall(newWall);
218  m_mesh->UpdateNodeOwningWalls(newWall);
219  }
220  else {
221  // This case requires the possible splitting of walls.
222  assert(m_walls.size() >= 2 && "There can't be one wall in the mesh.");
223 
224  auto& boundaryWalls = m_mesh->GetBoundaryPolygon()->GetWalls();
225 
226  if (m_mesh->GetNodeOwningNeighbors(node1).size() == 2) {
227  auto endpoint = oldBoundaryNodes.rbegin();
228  while (m_mesh->GetNodeOwningNeighbors(*endpoint).size() == 2) {
229  endpoint++;
230  assert(endpoint != oldBoundaryNodes.rend()
231  && "There should be at least two endpoints with more than 2 owners in the mesh.");
232  }
233 
234  auto wall = find_if(boundaryWalls.begin(), boundaryWalls.end(), [endpoint](Wall* w){return w->GetN1() == *endpoint;});
235  assert(wall != boundaryWalls.end() && "There should be a wall for this endpoint.");
236 
237  Wall* newWall = m_mesh->BuildWall(node1, (*wall)->GetN2(), (*wall)->GetC1(), (*wall)->GetC2());
238  newWall->WallAttributes::operator=(**wall);
239  newWall->GetC1()->AddWall(newWall);
240  newWall->GetC2()->AddWall(newWall);
241 
242  (*wall)->SetN2(node1);
243  (*wall)->SetLengthDirty();
244 
245  m_mesh->UpdateNodeOwningWalls(newWall);
246  }
247 
248  if (m_mesh->GetNodeOwningNeighbors(node2).size() == 2) {
249  auto endpoint = find(oldBoundaryNodes.begin(), oldBoundaryNodes.end(), node2);
250  while (m_mesh->GetNodeOwningNeighbors(*endpoint).size() == 2) {
251  endpoint++;
252  assert(endpoint != oldBoundaryNodes.end() && "There should be at least two endpoints with more than 2 owners in the mesh.");
253  }
254 
255  auto wall = find_if(boundaryWalls.begin(), boundaryWalls.end(), [endpoint](Wall* w){return w->GetN2() == *endpoint;});
256  assert(wall != boundaryWalls.end() && "There should be a wall for this endpoint.");
257 
258  Wall* newWall = m_mesh->BuildWall((*wall)->GetN1(), node2, (*wall)->GetC1(), (*wall)->GetC2());
259  newWall->WallAttributes::operator=(**wall);
260  newWall->GetC1()->AddWall(newWall);
261  newWall->GetC2()->AddWall(newWall);
262 
263  (*wall)->SetN1(node2);
264  (*wall)->SetLengthDirty();
265 
266  m_mesh->UpdateNodeOwningWalls(newWall);
267  }
268  }
269 
270  // Make a vector of node IDs for the new cell (in the right order).
271  list<Node*> cellNodes(oldBoundaryNodes.begin(), next(find(oldBoundaryNodes.begin(), oldBoundaryNodes.end(), node2)));
272  cellNodes.reverse();
273  cellNodes.push_front(newNode);
274 
275  vector<unsigned int> nodeIDs(cellNodes.size());
276  transform(cellNodes.begin(), cellNodes.end(), nodeIDs.begin(), [](Node* n){return n->GetIndex();});
277 
278  // Build the new cell.
279  auto newCell = m_mesh->BuildCell(nodeIDs);
280 
281  // Reassign the correct walls from the boundary polygon to the new cell.
282  const auto boundaryWalls = m_mesh->GetBoundaryPolygon()->GetWalls();
283  for (Wall* w : boundaryWalls) {
284  auto n1 = find(cellNodes.begin(), cellNodes.end(), w->GetN1());
285  auto n2 = find(cellNodes.begin(), cellNodes.end(), w->GetN2());
286  if (n1 != cellNodes.end() && n2 != cellNodes.end() && (*n2 != node1 || *n1 != node2)) {
287  w->SetC2(newCell);
288  m_mesh->GetBoundaryPolygon()->ReassignWall(w, newCell);
289  }
290  }
291 
292  // Create the new wall between the new cell and the boundary polygon.
293  Wall* newWall = m_mesh->BuildWall(node1, node2,
294  newCell, m_mesh->GetBoundaryPolygon());
295  newWall->GetC1()->AddWall(newWall);
296  newWall->GetC2()->AddWall(newWall);
297  m_mesh->UpdateNodeOwningWalls(newWall);
298 
299  // Update the boundary neighbors of the nodes in the new cell to have the correct neighbors.
300  for (Node* n : newCell->GetNodes()) {
301  if (n != newNode) {
302  auto& owners = m_mesh->GetNodeOwningNeighbors(n);
303  auto owner = find_if(owners.begin(), owners.end(),
304  [](const NeighborNodes& nb){return nb.GetCell()->GetIndex() == -1;});
305  assert(owner != owners.end()
306  && "Node without a boundary neighbor in new cell that wasn't new node.");
307  owners.erase(owner);
308  }
309  }
310 
311  auto boundary_cit = make_const_circular(newBoundaryNodes);
312  do {
313  m_mesh->GetNodeOwningNeighbors(*boundary_cit).push_back(
314  NeighborNodes(m_mesh->GetBoundaryPolygon(),
315  *prev(boundary_cit), *next(boundary_cit)));
316  } while (*boundary_cit++ != node2);
317 
318  //Construct the neighbor lists of the new cell.
319  m_mesh->ConstructNeighborList(newCell);
320  for (Cell* c : m_mesh->GetNeighbors(newCell)) {
321  m_mesh->ConstructNeighborList(c);
322  }
323  m_mesh->ConstructNeighborList(m_mesh->GetBoundaryPolygon());
324 
325  // Set the nodes of the cell to the boundary or not when necessary.
326  for (Node* n : newCell->GetNodes()) {
327  SetAtBoundary(n);
328  }
329 
330  assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
331  return newCell;
332  }
333 }
334 
335 void EditableMesh::DeleteCells(const list<Cell*>& cells)
336 {
337  assert(cells.size() < m_cells.size() && "Can't delete all cells in the mesh.");
338  assert(find_if(cells.begin(), cells.end(), [this](Cell* c) {return (find(m_cells.begin(), m_cells.end(), c) == m_cells.end());}) == cells.end()
339  && "There's a given cell that doesn't exist in the mesh.");
340 
341  if (cells.empty()) {
342  return;
343  }
344 
345  // Delete the given cells and update the walls.
346  auto cell = m_cells.begin();
347  while (cell != m_cells.end()) {
348  if (find(cells.begin(), cells.end(), *cell) != cells.end()) {
349  for (Node* n : (*cell)->GetNodes()) {
350  n->SetAtBoundary(true);
351  }
352 
353  for (Wall* w : (*cell)->GetWalls()) {
354  if (w->IsAtBoundary()) {
355  w->SetN1(nullptr);
356  w->SetN2(nullptr);
357  m_walls.erase(remove_if(m_walls.begin(), m_walls.end(),
358  [w](Wall* wall){return wall == w;}),
359  m_walls.end());
360  //m_walls.remove_if([w](Wall* wall){return wall == w;});
361  }
362  else if (w->GetC1() == *cell) {
363  w->SetC1(w->GetC2());
364  w->SetC2(m_mesh->GetBoundaryPolygon());
365 
366  Node* tmpNode = w->GetN1();
367  w->SetN1(w->GetN2());
368  w->SetN2(tmpNode);
369 
370  // No need to set the wall length dirty bit, as this is the same wall:
371  // The cells and nodes are just swapped so that C1 isn't the boundary polygon
372  }
373  else { //walls->GetC2() == cell
374  w->SetC2(m_mesh->GetBoundaryPolygon());
375  }
376  }
377 
378  m_mesh->GetCells().erase(cell);
379  }
380  else {
381  cell++;
382  }
383  }
384 
385  //If there is only one cell left, all walls need to be deleted.
386  if (m_cells.size() == 1) {
387  m_walls.clear();
388  m_cells.front()->GetWalls().clear();
389  m_mesh->GetBoundaryPolygon()->GetWalls().clear();
390  }
391 
392  for (auto n : m_nodes) {
393  m_mesh->GetNodeOwningNeighbors(n).clear();
394  }
395 
396  list<Node*> nodesToKeep;
397  for (auto c : m_cells) {
398  auto cit = make_const_circular(c->GetNodes(), c->GetNodes().begin());
399  do {
400  nodesToKeep.push_back(*cit);
401  m_mesh->GetNodeOwningNeighbors(*cit).push_back(
402  NeighborNodes(c, *prev(cit), *next(cit)));
403  } while (++cit != c->GetNodes().begin());
404  }
405 
406  // Delete the useless nodes.
407  list<Node*> boundaryNodes;
408  auto node = m_nodes.begin();
409  while (node != m_nodes.end()) {
410  if (find(nodesToKeep.begin(), nodesToKeep.end(), *node) == nodesToKeep.end()) {
411  m_mesh->RemoveNodeFromOwnerLists(*node);
412  m_nodes.erase(node);
413  }
414  else {
415  if ((*node)->IsAtBoundary()) {
416  boundaryNodes.push_back(*node);
417  }
418  node++;
419  }
420  }
421 
422  //Re-add the new boundaryNodes to the boundary polygon.
423  Node* current = boundaryNodes.front();
424  boundaryNodes.pop_front();
425  m_mesh->GetBoundaryPolygon()->GetNodes().clear();
426  m_mesh->GetBoundaryPolygon()->GetNodes().push_back(current);
427  while (!boundaryNodes.empty()) {
428  // These are to check whether there are two neighboring
429  // cells (which implies we have the wrong boundary node).
430  list<Node*> firstNeighbors;
431  // These are the real options, as the order of the nodes
432  // should be the same for all the cells (counterclockwise).
433  list<Node*> secondNeighbors;
434  for (NeighborNodes nb : m_mesh->GetNodeOwningNeighbors(current)) {
435  if (nb.GetNb1()->IsAtBoundary()) {
436  firstNeighbors.push_back(nb.GetNb1());
437  }
438  if (nb.GetNb2()->IsAtBoundary()) {
439  secondNeighbors.push_back(nb.GetNb2());
440  }
441  }
442 
443  auto next = find_if(boundaryNodes.begin(), boundaryNodes.end(),
444  [firstNeighbors, secondNeighbors](Node* n) {
445  return find(firstNeighbors.begin(), firstNeighbors.end(), n)
446  == firstNeighbors.end()
447  && find(secondNeighbors.begin(), secondNeighbors.end(), n)
448  != secondNeighbors.end();
449  });
450  assert(next != boundaryNodes.end()
451  && "Inconsistency: No boundary neighbor to current boundary node.");
452 
453  current = *next;
454  boundaryNodes.erase(next);
455  m_mesh->GetBoundaryPolygon()->GetNodes().push_back(current);
456  }
457 
458  auto cit = make_const_circular(m_mesh->GetBoundaryPolygon()->GetNodes(),
459  m_mesh->GetBoundaryPolygon()->GetNodes().begin());
460  do {
461  m_mesh->GetNodeOwningNeighbors(*cit).push_back(
462  NeighborNodes(m_mesh->GetBoundaryPolygon(), *prev(cit), *next(cit)));
463  } while (++cit != m_mesh->GetBoundaryPolygon()->GetNodes().begin());
464 
465  // Boundary walls.
466  auto cmp_boundary_walls = [](Wall* w1, Wall* w2) {
467  Cell* c1 = w1->GetC1()->GetIndex() != -1 ? w1->GetC1() : w1->GetC2();
468  Cell* c2 = w2->GetC1()->GetIndex() != -1 ? w2->GetC1() : w2->GetC2();
469  return c1->GetIndex() < c2->GetIndex();
470  };
471  auto eql_boundary_walls = [](Wall* w1, Wall* w2) {
472  Cell* c1 = w1->GetC1()->GetIndex() != -1 ? w1->GetC1() : w1->GetC2();
473  Cell* c2 = w2->GetC1()->GetIndex() != -1 ? w2->GetC1() : w2->GetC2();
474  return c1->GetIndex() == c2->GetIndex();
475  };
476 
477  list<Wall*> boundaryWalls;
478  for (auto w : m_walls) {
479  if (w->IsAtBoundary()) {
480  boundaryWalls.push_back(w);
481  }
482  }
483 
484  //TODO: Random merge of walls. Something better?
485  m_mesh->GetBoundaryPolygon()->GetWalls().clear();
486  boundaryWalls.sort(cmp_boundary_walls);
487  auto wall = adjacent_find(boundaryWalls.begin(), boundaryWalls.end(), eql_boundary_walls);
488  while (wall != boundaryWalls.end()) {
489  auto next = wall;
490  next++;
491  while (next != boundaryWalls.end() && eql_boundary_walls(*wall, *next)) {
492  bool merge = true;
493  if ((*wall)->GetN1() == (*next)->GetN1()) {
494  (*wall)->SetN1((*next)->GetN2());
495  }
496  else if ((*wall)->GetN1() == (*next)->GetN2()) {
497  (*wall)->SetN1((*next)->GetN1());
498  }
499  else if ((*wall)->GetN2() == (*next)->GetN1()) {
500  (*wall)->SetN2((*next)->GetN2());
501  }
502  else if ((*wall)->GetN2() == (*next)->GetN2()) {
503  (*wall)->SetN2((*next)->GetN1());
504  }
505  else {
506  merge = false;
507  }
508 
509  if (merge) {
510  (*wall)->SetLengthDirty();
511 
512  (*next)->GetC1()->GetWalls().remove(*next);
513  (*next)->GetC2()->GetWalls().remove(*next);
514 
515  Wall* tmpWall = *next;
516  boundaryWalls.erase(next);
517 
518  //Remove from the node owning walls list.
519  tmpWall->SetN1(nullptr);
520  tmpWall->SetN2(nullptr);
521  m_mesh->UpdateNodeOwningWalls(*wall);
522 
523  //m_walls.remove_if([tmpWall](Wall* w){return w == tmpWall;});
524  m_walls.erase(remove_if(
525  m_walls.begin(), m_walls.end(),
526  [tmpWall](Wall* w){return w == tmpWall;}),
527  m_walls.end());
528  next = wall;
529  }
530  next++;
531  }
532  wall = adjacent_find(++wall, boundaryWalls.end(), eql_boundary_walls);
533  }
534 
535  for (auto w : boundaryWalls) {
536  m_mesh->GetBoundaryPolygon()->GetWalls().push_back(w);
537  }
538 
539  m_mesh->GetBoundaryPolygon()->SetGeoDirty();
540 
541  m_mesh->ConstructNeighborList(m_mesh->GetBoundaryPolygon());
542 
543  for (auto c : m_mesh->GetCells()) {
544  m_mesh->ConstructNeighborList(c);
545  }
546 
547  // Fix the IDs.
548  FixNodeIDs();
549  FixWallIDs();
550  FixCellIDs();
551 
552  if (m_mesh->GetWalls().size() == 0) {
553  for (auto n : m_mesh->GetNodes()) {
554  m_mesh->GetNodeOwningWalls(n).clear();
555  }
556  }
557 
558  for (auto w : m_mesh->GetWalls()) {
559  m_mesh->UpdateNodeOwningWalls(w);
560  }
561 
562  assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
563 }
564 
566 {
567  list<QLineF> oldEdges;
568  list<QLineF> newEdges;
569  for (NeighborNodes neighbor : m_mesh->GetNodeOwningNeighbors(node)) {
570  QLineF oldEdge1((*neighbor.GetNb1())[0], (*neighbor.GetNb1())[1], (*node)[0], (*node)[1]);
571  QLineF oldEdge2((*neighbor.GetNb2())[0], (*neighbor.GetNb2())[1], (*node)[0], (*node)[1]);
572  QLineF newEdge((*neighbor.GetNb1())[0], (*neighbor.GetNb1())[1], (*neighbor.GetNb2())[0], (*neighbor.GetNb2())[1]);
573 
574  if (find(oldEdges.begin(), oldEdges.end(), oldEdge1) == oldEdges.end()) {
575  oldEdges.push_front(oldEdge1);
576  }
577  if (find(oldEdges.begin(), oldEdges.end(), oldEdge2) == oldEdges.end()) {
578  oldEdges.push_front(oldEdge2);
579  }
580  if (find(newEdges.begin(), newEdges.end(), newEdge) == newEdges.end()) {
581  newEdges.push_front(newEdge);
582  }
583  }
584 
585  if (IsDeletableNode(node) && IsConsistent(oldEdges, newEdges)){
586  DeleteTwoDegreeNodeInconsistent(node);
587  FixNodeIDs();
588 
589  assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
590 
591  return true;
592  }
593  else {
594  return false;
595  }
596 }
597 
598 bool EditableMesh::DeleteTwoDegreeNodeInconsistent(Node* node)
599 {
600  //Update owners of node
601  list<NeighborNodes> neighbors = m_mesh->GetNodeOwningNeighbors(node);
602 
603  if (neighbors.size() != 2) {
604  return false; //Tried to delete a node that isn't connected to exactly two other nodes.
605  }
606 
607  auto nodeOwningWalls = m_mesh->GetNodeOwningWalls(node);
608  if (nodeOwningWalls.size() == 1) {
609  nodeOwningWalls.front()->SetLengthDirty();
610  }
611  else {
612  assert(nodeOwningWalls.size() == 0 && m_walls.size() == 0
613  && "A two degree node with no wall owners is only possible for one cell in the mesh (no walls).");
614  }
615 
616  Node* n1 = neighbors.front().GetNb1();
617  Node* n2 = neighbors.front().GetNb2();
618  NeighborNodes nb1 = neighbors.front();
619  NeighborNodes nb2 = neighbors.back();
620 
621  nb1.GetCell()->GetNodes().erase(std::find(nb1.GetCell()->GetNodes().begin(), nb1.GetCell()->GetNodes().end(), node));
622  nb1.GetCell()->SetGeoDirty();
623  nb2.GetCell()->GetNodes().erase(std::find(nb2.GetCell()->GetNodes().begin(), nb2.GetCell()->GetNodes().end(), node));
624  nb2.GetCell()->SetGeoDirty();
625 
626  //Lambda to test whether the given neighbor contains 'node'.
627  auto is_nb = [node](const NeighborNodes& nb) {
628  return (nb.GetNb1() == node || nb.GetNb2() == node);
629  };
630 
631  //Update owners of the first neighbor of node.
632  {
633  list<NeighborNodes>& owners = m_mesh->GetNodeOwningNeighbors(n1);
634  auto owner = find_if(owners.begin(), owners.end(), is_nb);
635  while (owner != owners.end()) {
636  if (owner->GetNb1() == node) {
637  const NeighborNodes tmp(owner->GetCell(), n2, owner->GetNb2());
638  *owner = tmp;
639  }
640  else if (owner->GetNb2() == node) {
641  const NeighborNodes tmp(owner->GetCell(), owner->GetNb1(), n2);
642  *owner = tmp;
643  }
644  owner = find_if(owner, owners.end(), is_nb);
645  }
646  }
647 
648  //Update owners of the second neighbor of node.
649  {
650  list<NeighborNodes>& owners = m_mesh->GetNodeOwningNeighbors(n2);
651  auto owner = find_if(owners.begin(), owners.end(), is_nb);
652  while (owner != owners.end()) {
653  if (owner->GetNb1() == node) {
654  const NeighborNodes tmp(owner->GetCell(), n1, owner->GetNb2());
655  *owner = tmp;
656  }
657  else if (owner->GetNb2() == node) {
658  const NeighborNodes tmp(owner->GetCell(), owner->GetNb1(), n1);
659  *owner = tmp;
660  }
661  owner = find_if(owner, owners.end(), is_nb);
662  }
663  }
664 
665  m_mesh->RemoveNodeFromOwnerLists(node);
666  auto it = m_nodes.begin();
667  while (it != m_nodes.end()) {
668  if (node == *it) {
669  break;
670  }
671  it++;
672  }
673  m_nodes.erase(it);
674 
675  return true;
676 }
677 
678 bool EditableMesh::DisplaceNode(Node* node, double x, double y)
679 {
680  double old_x = (*node)[0];
681  double old_y = (*node)[1];
682  (*node)[0] = x;
683  (*node)[1] = y;
684 
685  if (IsConsistent(node)) {
686  for (const auto &nb : m_mesh->GetNodeOwningNeighbors(node)) {
687  nb.GetCell()->SetGeoDirty();
688  }
689  for (const auto &wall : m_mesh->GetNodeOwningWalls(node)) {
690  wall->SetLengthDirty();
691  }
692 
693  return true;
694  }
695  else {
696  (*node)[0] = old_x;
697  (*node)[1] = old_y;
698  return false;
699  }
700 }
701 
702 void EditableMesh::FixCellIDs()
703 {
704  auto cmp = [](Cell* c1, Cell* c2) {
705  return c1->GetIndex() < c2->GetIndex();
706  };
707 
708  sort(m_cells.begin(), m_cells.end(), cmp);
709  for (unsigned int i = 0; i < m_cells.size(); i++) {
710  m_cells[i]->m_index = i;
711  }
712 
713  m_mesh->m_cell_next_id = m_cells.size();
714 }
715 
716 void EditableMesh::FixNodeIDs()
717 {
718  auto cmp = [](Node* n1, Node* n2) {
719  return n1->GetIndex() < n2->GetIndex();
720  };
721 
722  sort(m_nodes.begin(), m_nodes.end(), cmp);
723  for (unsigned int i = 0; i < m_nodes.size(); i++) {
724  m_nodes[i]->m_index = i;
725  }
726 
727  m_mesh->m_node_next_id = m_nodes.size();
728 }
729 
730 void EditableMesh::FixWallIDs()
731 {
732  auto cmp = [](Wall* w1, Wall* w2) {
733  return w1->GetIndex() < w2->GetIndex();
734  };
735 
736  sort(m_walls.begin(), m_walls.end(),cmp);
737  for (unsigned int i = 0; i < m_walls.size(); i++) {
738  m_walls[i]->m_wall_index = i;
739  }
740 
741  m_mesh->m_wall_next_id = m_walls.size();
742 }
743 
744 std::shared_ptr<Mesh>& EditableMesh::GetMesh()
745 {
746  return m_mesh;
747 }
748 
749 bool EditableMesh::IsConsistent(const list<QLineF>& excluded, const list<QLineF>& included) const
750 {
751  // Lambda to check whether an edge is in the list of excluded edges.
752  auto is_excluded = [excluded](const QLineF& edge) {
753  return (find(excluded.begin(), excluded.end(), edge) != excluded.end());
754  };
755 
756  // Lambda to check whether an edge is in the list of included edges.
757  auto is_included = [included](const QLineF& edge) {
758  return (find(included.begin(), included.end(), edge) != included.end());
759  };
760 
761  // Include all the edges in the mesh.
762  list<QLineF> edges;
763  for (Cell* cell : m_cells) {
764  auto cit = make_const_circular(cell->GetNodes());
765  do {
766  QLineF newEdge((**cit)[0], (**cit)[1], (**next(cit))[0], (**next(cit))[1]);
767  if (!is_excluded(newEdge) && !is_included(newEdge)
768  && find(edges.begin(), edges.end(), newEdge) == edges.end()) {
769  edges.push_back(newEdge);
770  }
771  } while (++cit != cell->GetNodes().begin());
772  }
773 
774  // Check whether any of the included edges intersects with one
775  // of the mesh edges, excluding the edges in 'excluded'.
776  for (QLineF e1 : edges) {
777  for (QLineF e2 : included) {
778  QLineF::IntersectType intersection = e1.intersect(e2, nullptr);
779  if (intersection == QLineF::BoundedIntersection
780  && e1.p1() != e2.p1() && e1.p1() != e2.p2()
781  && e1.p2() != e2.p1() && e1.p2() != e2.p2()) {
782 
783  return false;
784  }
785  }
786  }
787 
788  return true;
789 }
790 
791 bool EditableMesh::IsConsistent(Node* node) const
792 {
793  list<Cell*> neighbors;
794  for (const NeighborNodes& nb : m_mesh->GetNodeOwningNeighbors(node)) {
795  if (nb.GetCell()->GetIndex() != -1 && find(neighbors.begin(), neighbors.end(), nb.GetCell()) == neighbors.end()) {
796  neighbors.push_back(nb.GetCell());
797  }
798  }
799 
800  if (node->IsAtBoundary()) {
801  for (Node* n : m_mesh->GetBoundaryPolygon()->GetNodes()) {
802  for (const NeighborNodes& nb : m_mesh->GetNodeOwningNeighbors(n)) {
803  if (nb.GetCell()->GetIndex() != -1 && find(neighbors.begin(), neighbors.end(), nb.GetCell()) == neighbors.end()) {
804  neighbors.push_back(nb.GetCell());
805  }
806  }
807  }
808  }
809 
810  list<QPolygonF> polygons;
811  for (Cell* c : neighbors) {
812  QPolygonF polygon;
813  for (Node* n : c->GetNodes()) {
814  polygon.append(QPointF((*n)[0], (*n)[1]));
815  }
816  polygons.push_back(polygon);
817  }
818 
819  for (auto p1 = polygons.begin(); p1 != polygons.end(); p1++) {
820  if (!PolygonUtils::IsSimplePolygon(*p1)) {
821  return false;
822  }
823 
824  for (auto p2 = next(p1); p2 != polygons.end(); p2++) {
825  if (p1->intersected(*p2).size() >= 3) {
826  return false;
827  }
828  }
829  }
830 
831  return true;
832 }
833 
835 {
836  if (!cell->HasBoundaryWall() || m_cells.size() <= 1
837  || find_if(m_cells.begin(), m_cells.end(), [cell](Cell* c){return c == cell;}) == m_cells.end()) {
838 
839  return false;
840  }
841 
842  for (Wall* wall : cell->GetWalls()) {
843  if (wall->IsAtBoundary()) {
844  auto nodes = cell->GetNodes();
845 
846  auto firstNode = (wall->GetC1() == cell ? wall->GetN2() : wall->GetN1());
847  rotate(nodes.begin(), next(find(nodes.begin(), nodes.end(), firstNode)), nodes.end());
848  auto secondNode = (wall->GetC1() == cell ? wall->GetN1() : wall->GetN2());
849 
850  for (auto it = nodes.begin(); it != find(nodes.begin(), nodes.end(), secondNode); it++) {
851  if ((*it)->IsAtBoundary()) {
852  //There's a node not on the wall that's at the boundary.
853  return false;
854  }
855  }
856 
857  return true;
858  }
859  }
860  assert(false && "There's no boundary wall, while this already has been verified.");
861 
862  return false;
863 }
864 
866 {
867  const list<NeighborNodes> neighbors = m_mesh->GetNodeOwningNeighbors(node);
868  if (neighbors.size() != 2) {
869  return false; //Tried to delete a node that isn't connected to exactly two other nodes.
870  }
871 
872  NeighborNodes nb1 = neighbors.front();
873  NeighborNodes nb2 = neighbors.back();
874 
875  assert(((nb1.GetNb1() == nb2.GetNb1() && nb1.GetNb2() == nb2.GetNb2())
876  || (nb1.GetNb1() == nb2.GetNb2() && nb1.GetNb2() == nb2.GetNb1()))
877  && nb1.GetCell() != nb2.GetCell()
878  && "Inconsistency in the neighbor list of the given node.");
879 
880  if (nb1.GetCell()->GetNodes().size() <= 3 || nb2.GetCell()->GetNodes().size() <= 3) {
881  return false; //After deleting the node, it wouldn't be a valid cell anymore.
882  }
883 
884  return true;
885 }
886 
887 void EditableMesh::ReplaceCell(Cell* cell, list<QPolygonF> newCells)
888 {
889  assert(find_if(newCells.begin(), newCells.end(),
890  [](const QPolygonF& p){return p.isClosed();}) == newCells.end()
891  && "All the given polygons should be open.");
892  assert(find_if(newCells.begin(), newCells.end(),
893  [](const QPolygonF& p){return p.size() < 3;}) == newCells.end()
894  && "There's a given polygon which is invalid. (less than three corners)");
895  assert(!newCells.empty() && "There are no new cells given.");
896  assert(all_of(newCells.begin(), newCells.end(),
897  [](const QPolygonF& p){return !PolygonUtils::IsClockwise(p);})
898  && "The given polygons weren't counter-clockwise.");
899  assert(all_of(cell->GetNodes().begin(), cell->GetNodes().end(), [newCells](Node* n){
900  for (const QPolygonF& p : newCells) {
901  if (p.contains(QPointF((*n)[0], (*n)[1]))) {
902  return true;
903  }
904  }
905  return false;
906  }) && "Not all the points of the cell were present in the new cells.");
907 
908  if (newCells.size() == 1) {
909  // newCells contains the old cell.
910  return;
911  }
912 
913  // Lambda to compare two points.
914  auto cmp_points = [](const QPointF& p1, const QPointF& p2) {
915  return p1 != p2 && (p1.x() < p2.x() || (p1.x() == p2.x() && p1.y() < p2.y()));
916  };
917 
918  // Lambda to compare two edges.
919  auto cmp_edges = [](const Edge& e1, const Edge& e2) {
920  return e1.GetFirst()->GetIndex() < e2.GetFirst()->GetIndex()
921  || (e1.GetFirst()->GetIndex() == e2.GetFirst()->GetIndex()
922  && e1.GetSecond()->GetIndex() < e2.GetSecond()->GetIndex());
923  };
924 
925  // Lambda to check whether to points are equal with a small perturbation.
926  auto is_perturbed = [](const QPointF& p1, const QPointF& p2) {
927  return fabs(p1.x() - p2.x()) < g_accuracy && fabs(p1.y() - p2.y()) < g_accuracy;
928  };
929 
930  // Lambda to perturb a line.
931  auto perturb_line = [](QLineF line) {
932  line.setLength(line.length() + g_accuracy);
933  line.setPoints(line.p2(), line.p1());
934  line.setLength(line.length() + g_accuracy);
935  line.setPoints(line.p2(), line.p1());
936  return line;
937  };
938 
939  // Delete the old cell and keep track of its nodes and its edge owners.
940  list<Cell*> cells;
941  map<QPointF, Node*, decltype(cmp_points)> newNodes(cmp_points);
942  map<Edge, Cell*, decltype(cmp_edges)> edgeOwners(cmp_edges);
943 
944 
945  auto edge_cit = make_const_circular(cell->GetNodes());
946  do {
947  Edge edge(*edge_cit, *next(edge_cit));
948  edgeOwners[edge] = m_mesh->FindEdgeNeighbor(cell, edge);
949  } while (++edge_cit != cell->GetNodes().begin());
950 
951  list<QPointF> points;
952  for (Node* node : cell->GetNodes()) {
953  points.push_back(QPointF((*node)[0], (*node)[1]));
954  newNodes[QPointF((*node)[0], (*node)[1])] = node;
955  }
956 
957  for (const QPolygonF& polygon : newCells) {
958  vector<unsigned int> nodeIDs;
959  for (const QPointF& point : polygon) {
960  points.push_back(point);
961 
962  if (newNodes.find(point) == newNodes.end()) {
963  newNodes[point] = m_mesh->BuildNode({{point.x(), point.y(), 0}}, false);
964  }
965  nodeIDs.push_back(newNodes[point]->GetIndex());
966  }
967 
968  Cell* c = m_mesh->BuildCell(nodeIDs);
969  c->CellAttributes::operator=(*cell);
970  cells.push_back(c);
971  }
972 
973  std::list<Wall*> newWalls;
974  Node* endpoint = nullptr; //Only used for the special case where we have no walls a priori.
975  for (const QPolygonF& polygon : newCells) {
976  auto polygonStd = polygon.toStdVector();
977  auto polygon_cit = make_const_circular(polygonStd);
978  do {
979  auto cell_cit = make_circular(cell->GetNodes());
980  do {
981  QPointF cell_p1((**cell_cit)[0], (**cell_cit)[1]);
982  QPointF cell_p2((**next(cell_cit))[0], (**next(cell_cit))[1]);
983  QLineF cell_line = perturb_line(QLineF(cell_p1, cell_p2));
984 
985  QLineF polygon_line = perturb_line(QLineF(*polygon_cit, *next(polygon_cit)));
986 
987  QPointF intersection;
988  if (polygon_line.intersect(cell_line, &intersection) == QLineF::BoundedIntersection) {
989  if (is_perturbed(intersection, *polygon_cit)) {
990  intersection = *polygon_cit;
991  assert(newNodes.find(intersection) != newNodes.end()
992  && "All intersections should already have been added.");
993  }
994  else if (is_perturbed(intersection, *next(polygon_cit))) {
995  intersection = *next(polygon_cit);
996  assert(newNodes.find(intersection) != newNodes.end()
997  && "All intersections should already have been added.");
998  }
999  else { //intersection with perturbed, parallel line in the middle.
1000  continue;
1001  }
1002 
1003  if (intersection != cell_p1 && intersection != cell_p2) {
1004  Edge edge(*cell_cit, *next(cell_cit));
1005  Edge edge1(*cell_cit, newNodes[intersection]);
1006  Edge edge2(newNodes[intersection], *next(cell_cit));
1007  Cell* neighbor = edgeOwners[edge];
1008 
1009  if (m_mesh->GetWalls().size() == 0 && endpoint == nullptr) {
1010  endpoint = newNodes[intersection];
1011  }
1012  else if (m_mesh->GetWalls().size() == 0) {
1013  Wall* newWall = m_mesh->BuildWall(newNodes[intersection], endpoint, cell, m_mesh->GetBoundaryPolygon());
1014  cell->AddWall(newWall);
1015  neighbor->AddWall(newWall);
1016  newWalls.push_back(newWall);
1017  newWall = m_mesh->BuildWall(endpoint, newNodes[intersection], cell, m_mesh->GetBoundaryPolygon());
1018  cell->AddWall(newWall);
1019  neighbor->AddWall(newWall);
1020  newWalls.push_back(newWall);
1021  }
1022  else {
1023  auto wall = find_if(cell->GetWalls().begin(), cell->GetWalls().end(), [edge](Wall* w){return w->IncludesEdge(edge);});
1024  assert(wall != cell->GetWalls().end()
1025  && find_if(next(wall), cell->GetWalls().end(), [edge](Wall* w){return w->IncludesEdge(edge);}) == cell->GetWalls().end()
1026  && "There should be exactly one wall associated with this edge.");
1027 
1028  // We know in which order the endpoints should be, as the edges are defined in the same order as the nodes in the deleted cell.
1029  Wall* newWall;
1030  if ((*wall)->GetC1() == cell) {
1031  newWall = m_mesh->BuildWall(newNodes[intersection], (*wall)->GetN2(), cell, neighbor);
1032  }
1033  else { // (*wall)->GetC2() == cell
1034  newWall = m_mesh->BuildWall(newNodes[intersection], (*wall)->GetN2(), neighbor, cell);
1035  }
1036  (*wall)->SetN2(newNodes[intersection]);
1037  (*wall)->SetLengthDirty();
1038  cell->AddWall(newWall);
1039  neighbor->AddWall(newWall);
1040 
1041  newWalls.push_back(*wall);
1042  newWalls.push_back(newWall);
1043  }
1044 
1045  edgeOwners[edge1] = neighbor;
1046  edgeOwners[edge2] = neighbor;
1047  edgeOwners.erase(edge);
1048 
1049  if (neighbor->GetIndex() != -1) {
1050  m_mesh->AddNodeToCell(neighbor, newNodes[intersection], *next(cell_cit), *cell_cit);
1051  }
1052  else {
1053  m_mesh->AddNodeToCell(neighbor, newNodes[intersection], *cell_cit, *next(cell_cit));
1054  }
1055  m_mesh->AddNodeToCell(cell, newNodes[intersection], *cell_cit, *next(cell_cit));
1056 
1057  cell_cit = cell->GetNodes().begin();
1058  }
1059  else if (newNodes.find(*polygon_cit) != newNodes.end() && newNodes.find(*next(polygon_cit)) != newNodes.end()) {
1060  Edge polygon_edge(newNodes[*polygon_cit], newNodes[*next(polygon_cit)]);
1061 
1062  if (!cell->HasEdge(polygon_edge)) {
1063  //When an edge of the new polygon goes right through a node that already existed in the old cell (and the edge isn't part of the old cell).
1064 
1065  Edge edge(*cell_cit, *next(cell_cit));
1066  Cell* neighbor = edgeOwners[edge];
1067 
1068  if (m_mesh->GetWalls().size() == 0 && endpoint == nullptr) {
1069  endpoint = newNodes[intersection];
1070  }
1071  else if (m_mesh->GetWalls().size() == 0) {
1072  if (endpoint != newNodes[intersection]) {
1073  Wall* newWall = m_mesh->BuildWall(newNodes[intersection], endpoint, cell, m_mesh->GetBoundaryPolygon());
1074  cell->AddWall(newWall);
1075  neighbor->AddWall(newWall);
1076  newWalls.push_back(newWall);
1077  newWall = m_mesh->BuildWall(endpoint, newNodes[intersection], cell, m_mesh->GetBoundaryPolygon());
1078  cell->AddWall(newWall);
1079  neighbor->AddWall(newWall);
1080  newWalls.push_back(newWall);
1081  }
1082  }
1083  else {
1084  auto wall = find_if(cell->GetWalls().begin(), cell->GetWalls().end(), [edge](Wall* w){return w->IncludesEdge(edge);});
1085  assert(wall != cell->GetWalls().end()
1086  && find_if(next(wall), cell->GetWalls().end(), [edge](Wall* w){return w->IncludesEdge(edge);}) == cell->GetWalls().end()
1087  && "There should be exactly one wall associated with this edge.");
1088 
1089  if ((*wall)->GetN1() != newNodes[intersection] && (*wall)->GetN2() != newNodes[intersection]) {
1090  Wall* newWall;
1091  if ((*wall)->GetC1() == cell) {
1092  newWall = m_mesh->BuildWall(newNodes[intersection], (*wall)->GetN2(), cell, neighbor);
1093  }
1094  else { // (*wall)->GetC2() == cell
1095  newWall = m_mesh->BuildWall(newNodes[intersection], (*wall)->GetN2(), neighbor, cell);
1096  }
1097  (*wall)->SetN2(newNodes[intersection]);
1098  (*wall)->SetLengthDirty();
1099  cell->AddWall(newWall);
1100  neighbor->AddWall(newWall);
1101 
1102  newWalls.push_back(*wall);
1103  newWalls.push_back(newWall);
1104  }
1105  }
1106  }
1107  }
1108  }
1109  } while (++cell_cit != cell->GetNodes().begin());
1110  } while (++polygon_cit != polygonStd.begin());
1111  }
1112 
1113  // Update neighbors.
1114  for (const auto& n : newNodes) {
1115  m_mesh->GetNodeOwningNeighbors(n.second).remove_if(
1116  [cell](const NeighborNodes& nb){return nb.GetCell() == cell;});
1117  }
1118 
1119  // Update wall owners.
1120  for (const auto& w : newWalls) {
1121  m_mesh->UpdateNodeOwningWalls(w);
1122  }
1123 
1124  // Create new walls in between new cells or reassign the walls from the old cell to the correct cells.
1125  map<Edge, Cell*, decltype(cmp_edges)> neighbors(cmp_edges);
1126  for (Cell* c : cells) {
1127  auto c_cit = make_const_circular(c->GetNodes());
1128  do {
1129  Edge edge(*c_cit, *next(c_cit));
1130  Edge reversedEdge(*next(c_cit), *c_cit);
1131 
1132  auto wall = find_if(cell->GetWalls().begin(), cell->GetWalls().end(),
1133  [edge](Wall* w){return w->IncludesEdge(edge);});
1134  if (wall != cell->GetWalls().end()) {
1135  assert(find_if(next(wall), cell->GetWalls().end(),
1136  [edge](Wall* w){return w->IncludesEdge(edge);})
1137  == cell->GetWalls().end()
1138  && "There should be exactly one wall associated with this edge.");
1139 
1140  if ((*wall)->GetC1() == cell) {
1141  (*wall)->SetC1(c);
1142  }
1143  else {
1144  (*wall)->SetC2(c);
1145  }
1146  cell->ReassignWall(*wall, c);
1147  m_mesh->UpdateNodeOwningWalls(*wall);
1148  }
1149  else if (edgeOwners.find(edge) == edgeOwners.end()
1150  && edgeOwners.find(reversedEdge) == edgeOwners.end()) {
1151  if (neighbors.find(reversedEdge) == neighbors.end()
1152  && neighbors.find(edge) == neighbors.end()) {
1153  neighbors[edge] = c;
1154  }
1155  else {
1156  Edge owner = (neighbors.find(edge) != neighbors.end()
1157  ? edge : reversedEdge);
1158 
1159  Cell* neighbor = neighbors[owner];
1160  neighbors.erase(owner);
1161 
1162  Wall* newWall = m_mesh->BuildWall(owner.GetFirst(),
1163  owner.GetSecond(), neighbor, c);
1164 
1165  c->AddWall(newWall);
1166  neighbor->AddWall(newWall);
1167  m_mesh->UpdateNodeOwningWalls(newWall);
1168  }
1169  }
1170  } while (++c_cit != c->GetNodes().begin());
1171  }
1172 
1173  m_mesh->RemoveFromNeighborList(cell);
1174 
1175  // Set the new nodes at the boundary if they should be.
1176  for (const auto& n : newNodes) {
1177  SetAtBoundary(n.second);
1178  }
1179 
1180  // Construct the neighbor lists.
1181  for (Cell* c : cells) {
1182  m_mesh->ConstructNeighborList(c);
1183  }
1184 
1185  list<Cell*> neighborCells;
1186  for (const auto& owner : edgeOwners) {
1187  auto neighbor = find(neighborCells.begin(), neighborCells.end(), owner.second);
1188  if (neighbor == neighborCells.end()) {
1189  m_mesh->ConstructNeighborList(owner.second);
1190  neighborCells.push_front(owner.second);
1191  }
1192  }
1193 
1194  auto meshCell = find_if(m_mesh->GetCells().begin(), m_mesh->GetCells().end(),
1195  [cell](Cell* c){return c == cell;});
1196  m_mesh->GetCells().erase(meshCell);
1197  FixCellIDs();
1198 
1199  assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
1200 }
1201 
1202 void EditableMesh::SetAtBoundary(Node* node)
1203 {
1204  auto boundary = find_if(
1205  m_mesh->GetNodeOwningNeighbors(node).begin(),
1206  m_mesh->GetNodeOwningNeighbors(node).end(),
1207  [](const NeighborNodes& nb){return nb.GetCell()->GetIndex() == -1;});
1208  node->SetAtBoundary(boundary != m_mesh->GetNodeOwningNeighbors(node).end());
1209 }
1210 
1212 {
1213  assert(m_mesh->GetEdgeOwners(edge).size() > 0 && "Edge doesn't exist in the mesh.");
1214 
1215  Wall* wall = m_mesh->FindWall(edge);
1216  Node* node = m_mesh->BuildNode({{((*edge.GetFirst())[0] + (*edge.GetSecond())[0]) / 2, ((*edge.GetFirst())[1] + (*edge.GetSecond())[1]) / 2, 0}}, m_mesh->IsAtBoundary(&edge));
1217  for (NeighborNodes neighbor : m_mesh->GetEdgeOwners(edge)) {
1218  if (find(neighbor.GetCell()->GetNodes().begin(), neighbor.GetCell()->GetNodes().end(), node) == neighbor.GetCell()->GetNodes().end()) {
1219  m_mesh->AddNodeToCell(neighbor.GetCell(), node, edge.GetFirst(), edge.GetSecond());
1220  }
1221  }
1222 
1223  if (m_mesh->GetWalls().size() > 0) {
1224  m_mesh->UpdateNodeOwningWalls(wall);
1225  }
1226 
1227  assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
1228  return node;
1229 }
1230 
1231 vector<Cell*> EditableMesh::SplitCell(Node* node1, Node* node2)
1232 {
1233  QLineF cut((*node1)[0], (*node1)[1], (*node2)[0], (*node2)[1]);
1234 
1235  list<Cell*> owners;
1236  for (auto nb : m_mesh->GetNodeOwningNeighbors(node1)) {
1237  if (nb.GetCell()->GetIndex() != -1) {
1238  owners.push_back(nb.GetCell());
1239  }
1240  }
1241  for (auto nb : m_mesh->GetNodeOwningNeighbors(node2)) {
1242  if (nb.GetCell()->GetIndex() != -1) {
1243  owners.push_back(nb.GetCell());
1244  }
1245  }
1246  owners.sort([](Cell* c1, Cell* c2) {return c1 < c2;});
1247 
1248  auto eql = [](Cell* c1, Cell* c2) {
1249  return c1 == c2;
1250  };
1251 
1252  list<Cell*> commonOwners;
1253  auto owner = adjacent_find(owners.begin(), owners.end(), eql);
1254  while (owner != owners.end()) {
1255  commonOwners.push_back(*owner);
1256  owner = adjacent_find(++owner, owners.end(), eql);
1257  }
1258 
1259  if (commonOwners.empty()) {
1260  return vector<Cell*>();
1261  }
1262 
1263  list<Cell*> possibilities;
1264  for (Cell* c : commonOwners) {
1265  QPolygonF polygon;
1266  for (Node* n : c->GetNodes()) {
1267  polygon.push_back(QPointF((*n)[0], (*n)[1]));
1268  }
1269  //if (polygon.containsPoint(center, Qt::OddEvenFill)) {
1270  auto polygons = PolygonUtils::SlicePolygon(polygon, cut);
1271  if (polygons.size() == 2) {
1272  possibilities.push_back(c);
1273  }
1274  }
1275  assert(possibilities.size() == 1 && "There can only be one cell that can be split.");
1276 
1277  auto cell = find_if(m_mesh->GetCells().begin(), m_mesh->GetCells().end(),
1278  [possibilities](Cell* c){return c == possibilities.front();});
1279 
1280  return SplitCell(*cell, node1, node2);
1281 }
1282 
1283 vector<Cell*> EditableMesh::SplitCell(Cell* cell, Node* node1, Node* node2)
1284 {
1285  assert(cell != nullptr && node1 != nullptr
1286  && node2 != nullptr && "The input isn't properly initialized.");
1287 
1288  vector<Cell*> newCells;
1289 
1290  // Input verification.
1291  if (!CanSplitCell(node1, node2, cell)) {
1292  return newCells;
1293  }
1294  else {
1295  for (NeighborNodes nb : m_mesh->GetNodeOwningNeighbors(node1)) {
1296  if (nb.GetNb1() == node2 || nb.GetNb2() == node2) {
1297  newCells.push_back(cell);
1298  return newCells;
1299  }
1300  }
1301 
1302  // Check whether the line between node1 and node2 intersects with
1303  // one of the edges of the cell (except for the endpoints).
1304  // If this is the case, return the old cell.
1305  QLineF cut((*node1)[0], (*node1)[1], (*node2)[0], (*node2)[1]);
1306  auto cit = make_const_circular(cell->GetNodes());
1307  do {
1308  QPointF point1((**cit)[0], (**cit)[1]);
1309  QPointF point2((**next(cit))[0], (**next(cit))[1]);
1310  QPointF intersection;
1311  if (cut.intersect(QLineF(point1, point2), &intersection) == QLineF::BoundedIntersection
1312  && intersection != point1 && intersection != point2) {
1313  newCells.push_back(cell);
1314  return newCells;
1315  }
1316  } while (++cit != cell->GetNodes().begin());
1317  }
1318 
1319  CellDivider cellDivider(m_mesh.get());
1320  for (Cell* c : cellDivider.GeometricSplitCell(cell, node1, node2)) {
1321  auto newCell = find_if(m_mesh->GetCells().begin(), m_mesh->GetCells().end(),
1322  [c](Cell* c2){ return c2 == c; });
1323  assert(newCell != m_mesh->GetCells().end() && "The new cell can't be found in the mesh.");
1324  newCells.push_back(*newCell);
1325  }
1326 
1327  assert(MeshCheck(*m_mesh).CheckAll() && "Mesh inconsistent.");
1328  return newCells;
1329 }
1330 
1331 } // namespace
bool IsDeletableNode(Node *node) const
Checks whether the given node is deletable.
STL namespace.
Interface for CellDivider.
A cell contains walls and nodes.
Definition: Cell.h:48
Namespace for SimPT tissue editor package.
Definition: Cell.h:32
Node in cell wall.
Definition: Node.h:39
Cell * GetCell() const
Return the cell of this Neighbor pair.
Combo header for circular iterator.
Interface for EditableMesh.
bool CheckAll() const
Runs all of the of checks to verify consistency of the mesh.
Definition: MeshCheck.cpp:63
const std::list< Wall * > & GetWalls() const
Access the cell's walls.
Definition: Cell.h:88
Interface for NeighborNodes.
static bool IsClockwise(const QPolygonF &polygon)
Checks whether the vertices of a polygon are ordered clockwise or counterclockwise.
Cell * CreateCell(Node *node1, Node *node2, const QPointF &newPoint)
Creates a new cell with the two given nodes and a third new point.
static std::list< QPolygonF > SlicePolygon(const QPolygonF &polygon, QLineF cut)
Slice an open, simple and valid polygon in multiple polygons with a given line.
CircularIterator< typename T::iterator > make_circular(T &c)
Helper produces non-const circular iterator whose range corresponds to the begin and end iterators of...
Namespace for the core simulator.
bool DeleteTwoDegreeNode(Node *node)
Deletes the given node from the mesh (which should stay consistent).
An Edge connects two nodes and is ambidextrous.
Definition: Edge.h:31
Interface for Cell.
bool IsDeletableCell(Cell *cell) const
Checks whether the given cell is deletable (the mesh should stay consistent).
Namespace for container related classes.
Interface for WallAttributes.
Checks mesh concistency.
Definition: MeshCheck.h:29
Interface for MeshCheck.
Interface for PBMBBuilder.
Class for handling cell division.
Definition: CellDivider.h:48
Class that directs ptree based mesh building process.
Definition: PBMBuilder.h:33
Node * GetNb1() const
Return first node of this Neighbor pair.
bool CanSplitCell(Node *node1, Node *node2, Cell *cell=nullptr) const
Checks whether the two given nodes can split a (given) cell (excluding the boundary polygon)...
std::shared_ptr< SimPT_Sim::Mesh > & GetMesh()
Get the mesh.
Structure of neighboring nodes: two neighboring nodes from standpoint of a given cell with an orderin...
Definition: NeighborNodes.h:34
void DeleteCells(const std::list< Cell * > &cells)
Deletes the given cell from the mesh.
Node * GetFirst() const
Get the first node of the edge.
Definition: Edge.h:45
const std::vector< Node * > & GetNodes() const
Access the nodes of cell's polygon.
Definition: Cell.h:79
std::vector< Cell * > SplitCell(Node *node1, Node *node2)
Does the same as SplitCell, but finds the correct cell that should be split.
Interface for Node.
int GetIndex() const
Return the index.
Definition: Cell.h:76
Node * GetNb2() const
Return second node of this Neighbor pair.
Interface for Edge.
bool DisplaceNode(Node *node, double x, double y)
Place a node at the given coordinates.
CircularIterator< T > next(CircularIterator< T > i)
Helper yields the position the iterator would have if moved forward (in circular fashion) by 1 positi...
Interface for CellAttributes.
ConstCircularIterator< typename T::const_iterator > make_const_circular(const T &c)
Helper produces const circular iterator whose range corresponds to the begin and end iterators of a c...
static bool IsSimplePolygon(const QPolygonF &polygon)
Checks whether a polygon is simple.
Interface for PolygonUtils.
Node * SplitEdge(const Edge &edge)
Split a given edge at the center.
EditableMesh(const boost::property_tree::ptree &pt)
Constructor.
CircularIterator< T > prev(CircularIterator< T > i)
Helper yields the position the iterator would have if moved backward (in circular fashion) by 1 posit...
A cell wall, runs between cell corner points and consists of wall elements.
Definition: Wall.h:48
Node * GetSecond() const
Get the second node of the edge.
Definition: Edge.h:51
Interface for Wall.
void ReplaceCell(Cell *cell, std::list< QPolygonF > newCells)
Replace a given cell with a set of new cells.
Interface for Mesh.