VPTissue Reference Manual
MeshCheck.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 "MeshCheck.h"
21 
22 #include "Cell.h"
23 #include "Edge.h"
24 #include "Mesh.h"
25 #include "NeighborNodes.h"
26 #include "Node.h"
27 
30 #include "util/misc/log_debug.h"
31 
32 #include <boost/property_tree/ptree.hpp>
33 #include <algorithm>
34 #include <array>
35 #include <cmath>
36 #include <functional>
37 #include <iomanip>
38 #include <ios>
39 #include <iostream>
40 #include <iterator>
41 #include <list>
42 #include <memory>
43 #include <set>
44 #include <sstream>
45 #include <string>
46 #include <tuple>
47 #include <vector>
48 
49 namespace SimPT_Sim { class Wall; }
50 
51 using namespace boost::property_tree;
52 using namespace std;
53 using namespace SimPT_Sim;
54 using namespace SimPT_Sim::Container;
55 
56 namespace SimPT_Sim {
57 
58 MeshCheck::MeshCheck(const Mesh& mesh)
59  : m_mesh(mesh)
60 {
61 }
62 
63 bool MeshCheck::CheckAll() const
64 {
65  auto overall = true;
66  overall = overall && CheckAreas();
67  overall = overall && CheckAtBoundaryNodes();
68  overall = overall && CheckCellBoundaryWallNodes();
69  overall = overall && CheckCellBoundaryWalls();
70  overall = overall && CheckCellIdsSequence();
71  overall = overall && CheckCellIdsUnique();
72  overall = overall && CheckEdgeOwners();
73  overall = overall && CheckMutuallyNeighbors();
74  overall = overall && CheckNodeIdsSequence();
75  overall = overall && CheckNodeIdsUnique();
76  overall = overall && CheckNodeOwningNeighbors();
77  overall = overall && CheckNodeOwningWalls();
78  overall = overall && CheckWallIdsSequence();
79  overall = overall && CheckWallIdsUnique();
80  overall = overall && CheckWallNeighborsList();
81  return overall;
82 }
83 
85 {
86  double area = 0.0;
87  for (const auto& c: m_mesh.GetCells()) {
88  area += c->GetArea();
89  }
90  const double delta = abs(area - m_mesh.GetBoundaryPolygon()->GetArea());
91 
92  const auto status = (delta < 1.0e-7);
93  assert( status && "Sum of cell areas disagrees with boundary polygon area!");
94 
95  return status;
96 }
97 
99 {
100  auto status = true;
101  for (const auto& node : m_mesh.GetNodes()) {
102  if (!(node->IsAtBoundary() == m_mesh.IsInBoundaryPolygon(node))) {
103  status = false;
104  break;
105  }
106  }
107  assert( status && "Nodes and boundary_polgon disagree on at_boundary status!");
108  return status;
109 }
110 
112 {
113  auto status = true;
114  list<Cell*> bw_cells;
115 
116  // Select cells with a boundary wall.
117  for (auto cell : m_mesh.GetCells()) {
118  if (cell->IsWallNeighbor(m_mesh.GetBoundaryPolygon())) {
119  bw_cells.push_back(cell);
120  }
121  }
122 
123  // For each cell, isolate the boundary wall and check that
124  // the intervening nodes are really on the boundary polygon.
125  for (auto cell : bw_cells) {
126  auto is_bw = [this](Wall* w) {
127  return this->m_mesh.IsAtBoundary(w);
128  };
129  const auto& walls = cell->GetWalls();
130 
131  //TODO: What if there are multiple boundary walls?
132  const auto it = find_if(walls.begin(), walls.end(), is_bw);
133  if (it == walls.end()) {
134  status = false;
135  break;
136  }
137  const auto bw = *it;
138  const Node* start = bw->GetN1();
139  const Node* stop = bw->GetN2();
140  if (cell == bw->GetC2()) {
141  swap(start, stop);
142  }
143  auto& nodes = cell->GetNodes();
144  const auto start_it = find(nodes.begin(), nodes.end(), start);
145  auto c_it = make_circular(nodes.begin(), nodes.end(), start_it);
146  do {
147  if (!m_mesh.IsInBoundaryPolygon(*c_it)) {
148  status = false;
149  break;
150  };
151  } while(*(++c_it) != stop);
152  if (!status) {
153  break;
154  }
155  }
156 
157  assert( status && "Boundary walls have nodes not on boundary!");
158  return status;
159 }
160 
162 {
163  auto status = true;
164  if (m_mesh.GetCells().size() > 1U) {
165  for (auto cell : m_mesh.GetCells()) {
166  if (cell->IsWallNeighbor(m_mesh.GetBoundaryPolygon()) != cell->HasBoundaryWall()) {
167  status = false;
168  break;
169  }
170  }
171  }
172  assert( status && "Node and wall criteria disagree in at_boundary status!");
173  return status;
174 }
175 
177 {
178  bool status = true;
179 
180  auto cells = m_mesh.GetCells();
181  for(unsigned int i = 0; i < cells.size(); ++i) {
182  if (static_cast<unsigned int>(cells[i]->GetIndex()) != i) {
183  status = false;
184  break;
185  }
186  }
187 
188  assert( status && "Cell Ids not in sequence.");
189  return status;
190 }
191 
193 {
194  bool status = true;
195 
196  for (const auto& c1: m_mesh.GetCells()) {
197  for (const auto& c2: m_mesh.GetCells()) {
198  if (c1->GetIndex() == c2->GetIndex() && c1 != c2) {
199  status = false;
200  break;
201  }
202  }
203  if (status == false) {
204  break;
205  }
206  };
207 
208  assert( status && "Cell Ids not unique.");
209  return status;
210 }
212 {
213  auto status = true;
214  for (auto cell : m_mesh.GetCells()) {
215  const auto& nodes = cell->GetNodes();
216  const auto& walls = cell->GetWalls();
217  for (const auto& wall : walls) {
218  const auto cit1 = find(nodes.begin(), nodes.end(), wall->GetN1());
219  const auto cit2 = find(nodes.begin(), nodes.end(), wall->GetN2());
220  const auto cit_start = (wall->GetC1() == cell) ? cit1 : cit2;
221  const auto cit_stop = (wall->GetC1() == cell) ? cit2 : cit1;
222  const auto nghb_cell = (wall->GetC1() == cell) ? wall->GetC2() : wall->GetC1();
223 
224  auto cit = make_const_circular(nodes, cit_start);
225  do {
226  const Edge e(*cit, *next(cit));
227  const auto edge_owners = m_mesh.GetEdgeOwners(e);
228  status = cell->HasEdge(e) && nghb_cell->HasEdge(e);
229  status = status && (edge_owners.size() == 4 );
230  for (const auto& owner : edge_owners) {
231  status = status &&
232  (owner.GetCell() == cell || owner.GetCell() == nghb_cell);
233  }
234  if (!status) {
235  break;
236  }
237  } while(++cit != cit_stop);
238 
239  if (!status) {
240  break;
241  }
242  }
243  if (!status) {
244  break;
245  }
246  }
247 
248  assert( status && "EdgeOwners relation not OK!");
249  return status;
250 }
251 
253 {
254  auto status = true;
255  for (auto cell1 : m_mesh.GetCells()) {
256  for (auto cell2 : m_mesh.GetCells()) {
257  if (cell1->IsWallNeighbor(cell2) != cell2->IsWallNeighbor(cell1)) {
258  status = false;
259  break;
260  }
261  }
262  }
263 
264  assert( status && "Neighbor relation is not symmetric!");
265  return status;
266 }
267 
269 {
270  bool status = true;
271 
272  const auto& nodes = m_mesh.GetNodes();
273  for(unsigned int i = 0; i < nodes.size(); ++i) {
274  if (static_cast<unsigned int>(nodes[i]->GetIndex()) != i) {
275  status = false;
276  break;
277  }
278  }
279 
280  assert( status && "Cell Ids not in sequence.");
281  return status;
282 }
283 
285 {
286  bool status = true;
287  for (const auto& n1: m_mesh.GetNodes()) {
288  for (const auto& n2: m_mesh.GetNodes()) {
289  if (n1->GetIndex() == n2->GetIndex() && n1 != n2) {
290  status = false;
291  break;
292  }
293  }
294  if (status == false) {
295  break;
296  }
297  };
298 
299  assert( status && "Node Ids not unique.");
300  return status;
301 }
302 
304 {
305  auto status = true;
306  for (auto cell : m_mesh.GetCells()) {
307  const auto& nodes = cell->GetNodes();
308  auto cit_node = make_const_circular(nodes, nodes.begin());
309  do {
310  NeighborNodes nn(cell, *prev(cit_node), *next(cit_node));
311  const auto& owner = m_mesh.GetNodeOwningNeighbors(*cit_node);
312  const auto owner_it = find(owner.begin(), owner.end(), nn);
313  const auto owner_it2 = find(next(owner_it), owner.end(), nn);
314  status = status && (owner_it != owner.end()) && (owner_it2 == owner.end());
315  } while(++cit_node != nodes.begin());
316  }
317 
318  assert( status && "NodeOwningNeighbors relation inconsistent!");
319  return status;
320 }
321 
323 {
324  auto status = true;
325 
326  // Make a copy of the node owning walls so that we can remove every wall
327  // and check afterwards whether the map is empty
328  auto n_o_w_copy = m_mesh.m_node_owning_walls;
329  for (const auto& wall : m_mesh.GetWalls()) {
330  const auto& nodes = wall->GetC1()->GetNodes();
331  auto from = make_const_circular(nodes, find(nodes.begin(), nodes.end(), wall->GetN1()));
332  auto to = make_const_circular(nodes, find(nodes.begin(), nodes.end(), wall->GetN2()));
333 
334  // Check whether the endpoints of the wall are contained in C1
335  assert(from != nodes.end() && "Couldn't find N1 in C1, wall inconsistent");
336  assert(to != nodes.end() && "Couldn't find N2 in C1, wall inconsistent");
337  if (from == nodes.end() || to == nodes.end()) {
338  status = false;
339  continue;
340  }
341 
342  // Check whether this wall is in the node owning walls of all his nodes
343  do {
344  auto& owning_walls = n_o_w_copy[*from];
345  const auto it = std::find(owning_walls.begin(), owning_walls.end(), wall);
346  if (it == owning_walls.end()) {
347  status = false;
348  } else {
349  owning_walls.erase(it);
350  }
351  } while (from++ != to);
352  }
353 
354  assert(status && "Not every wall is associated with all its nodes");
355 
356  // Now check whether there are no more entries in the node owning walls than the ones we just found
357  for (const auto& pair : n_o_w_copy) {
358  if (!pair.second.empty()) {
359  status = false;
360  }
361  }
362 
363  assert(status && "Some node has an owning walls of which it is not part");
364  return status;
365 }
366 
368 {
369  bool status = true;
370 
371  const auto& walls = m_mesh.GetWalls();
372  for(unsigned int i = 0; i < walls.size(); ++i) {
373  if (static_cast<unsigned int>(walls[i]->GetIndex()) != i) {
374  status = false;
375  break;
376  }
377  }
378 
379  assert( status && "Wall Ids not in sequence.");
380  return status;
381 }
382 
384 {
385  bool status = true;
386  for (const auto& w1: m_mesh.GetNodes()) {
387  for (const auto& w2: m_mesh.GetNodes()) {
388  if (w1->GetIndex() == w2->GetIndex() && w1 != w2) {
389  status = false;
390  break;
391  }
392  }
393  if (status == false) {
394  break;
395  }
396  };
397 
398  assert( status && "Wall Ids not unique.");
399  return status;
400 }
401 
403 {
404  auto status = true;
405  for (auto cell1 : m_mesh.GetCells()) {
406  for (auto cell2 : m_mesh.GetCells()) {
407  if (cell1->IsWallNeighbor(cell2)) {
408  if (m_mesh.m_wall_neighbors.count(cell1) != 1) {
409  status = false;
410  break;
411  }
412  const auto& nb_list = m_mesh.m_wall_neighbors.at(cell1);
413  const auto it = find(nb_list.begin(), nb_list.end(), cell2);
414  if (it == nb_list.end()) {
415  status = false;
416  break;
417  }
418  }
419  }
420  }
421 
422  assert( status && "Wall_neighbor list in mesh has inconsistent info!");
423  return status;
424 }
425 
426 }
bool CheckNodeIdsSequence() const
Verifies that node ids correspond to storage sequence.
Definition: MeshCheck.cpp:268
bool IsAtBoundary(const Edge *e) const
True iff the edge is on the mesh boundary.
Definition: Mesh.cpp:636
STL namespace.
Node in cell wall.
Definition: Node.h:39
bool CheckCellBoundaryWalls() const
Verifies that for all cells the Cell::HasBoundaryWall agrees with Mesh::IsWallNeighbor for boundary p...
Definition: MeshCheck.cpp:161
bool CheckAll() const
Runs all of the of checks to verify consistency of the mesh.
Definition: MeshCheck.cpp:63
Interface for NeighborNodes.
bool CheckNodeIdsUnique() const
Verifies that node identifiers are unique.
Definition: MeshCheck.cpp:284
const std::vector< Cell * > & GetCells() const
The cells of this mesh, EXCLUDING the boundary polygon.
Definition: Mesh.h:108
bool CheckMutuallyNeighbors() const
Verifies that Cell::IsWallNeighbor is a symmetric relationship.
Definition: MeshCheck.cpp:252
bool CheckWallIdsUnique() const
Verifies that wall identifiers are unique.
Definition: MeshCheck.cpp:383
CircularIterator< typename T::iterator > make_circular(T &c)
Helper produces non-const circular iterator whose range corresponds to the begin and end iterators of...
bool CheckAtBoundaryNodes() const
Verifies Node::IsAtBoundary is true iff node is part of boundary polygon, for every node in the mesh...
Definition: MeshCheck.cpp:98
Namespace for the core simulator.
An Edge connects two nodes and is ambidextrous.
Definition: Edge.h:31
bool CheckCellIdsUnique() const
Verifies that cell identifiers are unique.
Definition: MeshCheck.cpp:192
Interface for Cell.
std::vector< NeighborNodes > GetEdgeOwners(const Edge &edge) const
Find the intersecting neighbors of the end points of the edge.
Definition: Mesh.cpp:455
Namespace for container related classes.
Interface for MeshCheck.
Macro defs for debug and logging.
bool CheckAreas() const
Verifies that sum of all cell areas equals the area of the boundary polygon.
Definition: MeshCheck.cpp:84
bool CheckWallNeighborsList() const
Verifies whether every cell has exactly one entry in the wall neighbor list (except when there is onl...
Definition: MeshCheck.cpp:402
Structure of neighboring nodes: two neighboring nodes from standpoint of a given cell with an orderin...
Definition: NeighborNodes.h:34
Cell * GetBoundaryPolygon() const
Get the boundary polygon of the mesh.
Definition: Mesh.h:105
bool IsInBoundaryPolygon(Node *n) const
True iff the node is one of the nodes of the boundary polygon.
Definition: Mesh.cpp:657
Interface for Node.
ConstCircularIterator class and helper functions.
const std::vector< Node * > & GetNodes() const
The nodes of the mesh.
Definition: Mesh.h:136
std::list< NeighborNodes > & GetNodeOwningNeighbors(Node *n)
Get the neighbors associated with this node.
Definition: Mesh.h:124
Interface for Edge.
Structure of cells; key data structure.
Definition: Mesh.h:62
CircularIterator< T > next(CircularIterator< T > i)
Helper yields the position the iterator would have if moved forward (in circular fashion) by 1 positi...
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...
const std::vector< Wall * > & GetWalls() const
The walls of this mesh.
Definition: Mesh.h:148
bool CheckCellBoundaryWallNodes() const
Verifies that for each cell, each boundary wall (other cell owning the wall is the boundary polygon) ...
Definition: MeshCheck.cpp:111
CircularIterator< T > prev(CircularIterator< T > i)
Helper yields the position the iterator would have if moved backward (in circular fashion) by 1 posit...
CircularIterator class and helper functions.
bool CheckCellIdsSequence() const
Verifies that cell ids correspond to storage sequence.
Definition: MeshCheck.cpp:176
bool CheckNodeOwningNeighbors() const
Verifies that for each node its NeighborNodes constructs appear once and only once in the NodeOwningN...
Definition: MeshCheck.cpp:303
bool CheckWallIdsSequence() const
Verifies that wall ids correspond to storage sequence.
Definition: MeshCheck.cpp:367
bool CheckNodeOwningWalls() const
Verifies that the owning walls of every node are exactly those walls to which the node belongs...
Definition: MeshCheck.cpp:322
bool CheckEdgeOwners() const
Verifies that every edge has exactly two neighbor cells by using Mesh::GetEdgeOwners.
Definition: MeshCheck.cpp:211
double GetArea() const
Return the area of the cell.
Definition: Cell.cpp:178
A cell wall, runs between cell corner points and consists of wall elements.
Definition: Wall.h:48
Interface for Mesh.