Chaste  Release::2018.1
VoronoiVertexMeshGenerator.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "VoronoiVertexMeshGenerator.hpp"
37 
38 #if BOOST_VERSION >= 105200
39 
40 using boost::polygon::voronoi_builder;
41 using boost::polygon::voronoi_diagram;
42 typedef boost::polygon::point_data<int> boost_point;
43 
45  unsigned numElementsY,
46  unsigned numRelaxationSteps,
47  double elementTargetArea)
48  : mpMesh(nullptr),
49  mpTorMesh(nullptr),
50  mNumElementsX(numElementsX),
51  mNumElementsY(numElementsY),
52  mNumRelaxationSteps(numRelaxationSteps),
53  mElementTargetArea(elementTargetArea),
54  mTol(1e-7),
55  mMaxExpectedNumSidesPerPolygon(17)
56 {
57  this->ValidateInputAndSetMembers();
58  this->GenerateVoronoiMesh();
59 }
60 
61 VoronoiVertexMeshGenerator::~VoronoiVertexMeshGenerator()
62 {
63  if (mpMesh)
64  {
65  delete mpMesh;
66  }
67  if (mpTorMesh)
68  {
69  delete mpTorMesh;
70  }
71 }
72 
73 void VoronoiVertexMeshGenerator::GenerateVoronoiMesh()
74 {
81  // Get initial seed locations
82  std::vector<c_vector<double, 2> > seed_locations = this->GetInitialPointLocations();
83  this->ValidateSeedLocations(seed_locations);
84 
85  // We now create the initial Voronoi tessellation. This method updates mpMesh.
86  this->CreateVoronoiTessellation(seed_locations);
87 
92  for (unsigned relaxation = 0 ; relaxation < mNumRelaxationSteps ; relaxation++)
93  {
94  seed_locations.clear();
95  seed_locations = this->GetElementCentroidsFromMesh();
96 
97  this->CreateVoronoiTessellation(seed_locations);
98  }
99 
100  // We need to modify the node locations to achieve the correct target average element area
101  double scale_factor = double(mMaxNumElems) * sqrt(mElementTargetArea);
102  for (unsigned node_idx = 0 ; node_idx < mpMesh->GetNumNodes() ; node_idx++)
103  {
104  c_vector<double, 2>& node_location = mpMesh->GetNode(node_idx)->rGetModifiableLocation();
105  node_location *= scale_factor;
106  }
107 }
108 
109 MutableVertexMesh<2,2>* VoronoiVertexMeshGenerator::GetMesh()
110 {
111  return mpMesh;
112 }
113 
114 MutableVertexMesh<2,2>* VoronoiVertexMeshGenerator::GetMeshAfterReMesh()
115 {
116  mpMesh->ReMesh();
117  return mpMesh;
118 }
119 
120 Toroidal2dVertexMesh* VoronoiVertexMeshGenerator::GetToroidalMesh()
121 {
122  /*
123  * METHOD OUTLINE:
124  *
125  * 1. Copy nodes and elements from mpMesh into a new MutableVertexMesh, so no data is shared between mpMesh and the
126  * nodes and elements we will be working with. There are no available copy constructors, so this is done from
127  * scratch.
128  *
129  * 2. Identify which nodes on the boundary are congruent to each other. This is done by recursively using the
130  * helper function CheckForCongruentNodes().
131  *
132  * 3. Create mpTorMesh by copying the subset of remaining nodes and all elements in the new MutableVertexMesh.
133  * Some elements have been modified to replace nodes by congruent partner nodes and some nodes have been deleted.
134  */
135 
136  // The width and height of the mesh for periodicity purposes
137  double width = mNumElementsX * sqrt(mElementTargetArea);
138  double height = mNumElementsY * sqrt(mElementTargetArea);
139 
140  // We need to construct new nodes and elements so we don't have mpTorMesh sharing data with mpMesh
141  std::vector<Node<2>*> new_nodes(mpMesh->GetNumNodes());
142  std::vector<VertexElement<2,2>*> new_elems(mpMesh->GetNumElements());
143 
144  // Copy nodes
145  for (unsigned node_counter = 0 ; node_counter < mpMesh->GetNumNodes() ; node_counter++)
146  {
147  Node<2>* p_node_to_copy = mpMesh->GetNode(node_counter);
148 
149  // Get all the information about the node we are copying
150  unsigned copy_index = p_node_to_copy->GetIndex();
151  bool copy_is_boundary = p_node_to_copy->IsBoundaryNode();
152  const c_vector<double, 2>& copy_location = p_node_to_copy->rGetLocation();
153 
154  // There should not be any 'gaps' in node numbering, but we will assert just to make sure
155  assert(copy_index < mpMesh->GetNumNodes());
156 
157  // Create a new node and place it in index order. Every node in a periodic mesh is non-boundary.
158  new_nodes[copy_index] = new Node<2>(copy_index, copy_location, copy_is_boundary);
159  }
160 
161  // Copy elements
162  for (unsigned elem_counter = 0; elem_counter < mpMesh->GetNumElements(); elem_counter++)
163  {
164  VertexElement<2,2>* p_elem_to_copy = mpMesh->GetElement(elem_counter);
165 
166  // Get the information relating to the element we are copying
167  unsigned copy_index = p_elem_to_copy->GetIndex();
168  unsigned copy_num_nodes = p_elem_to_copy->GetNumNodes();
169 
170  // There should not be any 'gaps' in element numbering, but we will assert just to make sure
171  assert(copy_index < mpMesh->GetNumElements());
172 
173  // The vertex element is created from a vector of nodes
174  std::vector<Node<2>*> nodes_this_elem;
175 
176  // Loop through the nodes in p_elem_to_copy and add the corresponding nodes that we have already copied
177  for (unsigned node_local_idx = 0 ; node_local_idx < copy_num_nodes ; node_local_idx++)
178  {
179  Node<2>* p_local_node = p_elem_to_copy->GetNode(node_local_idx);
180  unsigned local_node_global_idx = p_local_node->GetIndex();
181  nodes_this_elem.push_back(new_nodes[local_node_global_idx]);
182  }
183 
184  // Create a new node and place it in index order
185  new_elems[copy_index] = new VertexElement<2,2>(copy_index, nodes_this_elem);
186  }
187 
188  // We can now create the mesh with new_elements and the subset of new_nodes
189  MutableVertexMesh<2,2>* p_temp_mesh = new MutableVertexMesh<2,2>(new_nodes, new_elems);
190 
191  /*
192  * Recursively associate congruent nodes.
193  *
194  * For each node currently on the boundary, we loop through all other boundary nodes and check against all eight
195  * possible congruent locations. If any one coincides, we replace the identified node with its congruent partner,
196  * delete and remove the identified node from the mesh, and start again from scratch.
197  *
198  * If a node has no congruent partners, we simply mark it as non-boundary, as it is then a node that needs to appear
199  * in the final toroidal mesh.
200  *
201  * This recursion is guaranteed to terminate as the number of boundary nodes decreases by precisely one each time
202  * CheckForCongruentNodes() is called, and CheckForCongruentNodes() returns false if there are no boundary nodes.
203  */
204  bool re_check = true;
205 
206  while (re_check)
207  {
208  re_check = this->CheckForCongruentNodes(p_temp_mesh, width, height);
209  }
210 
211  // We now copy the nodes and elements into a new toroidal mesh, and delete p_temp_mesh
212  new_nodes.clear();
213  new_nodes.resize(p_temp_mesh->GetNumNodes());
214 
215  new_elems.clear();
216  new_elems.resize(p_temp_mesh->GetNumElements());
217 
218  // Copy nodes
219  for (unsigned node_counter = 0 ; node_counter < p_temp_mesh->GetNumNodes() ; node_counter++)
220  {
221  Node<2>* p_node_to_copy = p_temp_mesh->GetNode(node_counter);
222 
223  // Get all the information about the node we are copying
224  unsigned copy_index = p_node_to_copy->GetIndex();
225  const c_vector<double, 2>& copy_location = p_node_to_copy->rGetLocation();
226 
227  // No nodes should be boundary nodes
228  assert(!p_node_to_copy->IsBoundaryNode());
229 
230  // There should not be any 'gaps' in node numbering, but we will assert just to make sure
231  assert(copy_index < p_temp_mesh->GetNumNodes());
232 
233  // Create a new node and place it in index order. Every node in a periodic mesh is non-boundary.
234  new_nodes[copy_index] = new Node<2>(copy_index, copy_location, false);
235  }
236 
237  // Copy elements
238  for (unsigned elem_counter = 0; elem_counter < p_temp_mesh->GetNumElements(); elem_counter++)
239  {
240  VertexElement<2,2>* p_elem_to_copy = p_temp_mesh->GetElement(elem_counter);
241 
242  // Get the information relating to the element we are copying
243  unsigned copy_index = p_elem_to_copy->GetIndex();
244  unsigned copy_num_nodes = p_elem_to_copy->GetNumNodes();
245 
246  // There should not be any 'gaps' in element numbering, but we will assert just to make sure
247  assert(copy_index < p_temp_mesh->GetNumElements());
248 
249  // The vertex element is created from a vector of nodes
250  std::vector<Node<2>*> nodes_this_elem;
251 
252  // Loop through the nodes in p_elem_to_copy and add the corresponding nodes that we have already copied
253  for (unsigned node_local_idx = 0 ; node_local_idx < copy_num_nodes ; node_local_idx++)
254  {
255  Node<2>* p_local_node = p_elem_to_copy->GetNode(node_local_idx);
256 
257  unsigned local_node_global_idx = p_local_node->GetIndex();
258 
259  nodes_this_elem.push_back(new_nodes[local_node_global_idx]);
260  }
261 
262  // Create a new node and place it in index order
263  new_elems[copy_index] = new VertexElement<2,2>(copy_index, nodes_this_elem);
264  }
265 
266  delete p_temp_mesh;
267 
268  /*
269  * We now create the mesh with new_elements and new_nodes. We immediately call ReMesh() to tidy up any short edges.
270  *
271  * We then reposition all nodes to be within the box [0, width]x[0, height] to make mesh look nicer when visualised.
272  */
273  mpTorMesh = new Toroidal2dVertexMesh(width, height, new_nodes, new_elems);
274  mpTorMesh->ReMesh();
275 
276  c_vector<double, 2> min_x_y;
277  min_x_y[0] = DBL_MAX;
278  min_x_y[1] = DBL_MAX;
279 
280  // First loop is to calculate the correct offset, min_x_y
281  for (unsigned node_idx = 0 ; node_idx < mpTorMesh->GetNumNodes() ; node_idx++)
282  {
283  const c_vector<double, 2>& r_this_node_location = mpTorMesh->GetNode(node_idx)->rGetLocation();
284 
285  if (r_this_node_location[0] < min_x_y[0])
286  {
287  min_x_y[0] = r_this_node_location[0];
288  }
289  if (r_this_node_location[1] < min_x_y[1])
290  {
291  min_x_y[1] = r_this_node_location[1];
292  }
293  }
294 
295  // Second loop applies the offset, min_x_y, to each node in the mesh
296  for (unsigned node_idx = 0 ; node_idx < mpTorMesh->GetNumNodes() ; node_idx++)
297  {
298  mpTorMesh->GetNode(node_idx)->rGetModifiableLocation() -= min_x_y;
299  }
300 
301  // Third loop is to reposition any nodes that are now outside the bounding rectangle
302  for (unsigned node_idx = 0 ; node_idx < mpTorMesh->GetNumNodes() ; node_idx++)
303  {
304  if (mpTorMesh->GetNode(node_idx)->rGetLocation()[0] >= width)
305  {
306  mpTorMesh->GetNode(node_idx)->rGetModifiableLocation()[0] -= width;
307  }
308  if (mpTorMesh->GetNode(node_idx)->rGetLocation()[1] >= height)
309  {
310  mpTorMesh->GetNode(node_idx)->rGetModifiableLocation()[1] -= height;
311  }
312  }
313 
314  return mpTorMesh;
315 }
316 
317 bool VoronoiVertexMeshGenerator::CheckForCongruentNodes(MutableVertexMesh<2,2>* pMesh, double width, double height)
318 {
319  // First find all the current boundary nodes in pMesh
320  std::vector<Node<2>*> boundary_nodes;
321  for (unsigned node_idx = 0 ; node_idx < pMesh->GetNumNodes() ; node_idx++)
322  {
323  Node<2>* p_node = pMesh->GetNode(node_idx);
324 
325  if (p_node->IsBoundaryNode())
326  {
327  boundary_nodes.push_back(p_node);
328  }
329  }
330 
331  // If there are no boundary nodes, return false (we're done checking congruencies)
332  if (boundary_nodes.empty())
333  {
334  return false;
335  }
336 
337  // Otherwise, calculate the eight possible congruent locations for the current node
338  Node<2>* p_node_a = *(boundary_nodes.begin());
339  c_vector<double, 2> node_a_pos = p_node_a->rGetLocation();
340  std::vector<c_vector<double,2> > congruent_locations(8, node_a_pos);
341 
342  congruent_locations[0][0] += width;
343 
344  congruent_locations[1][1] += height;
345 
346  congruent_locations[2][0] -= width;
347 
348  congruent_locations[3][1] -= height;
349 
350  congruent_locations[4][0] += width;
351  congruent_locations[4][1] += height;
352 
353  congruent_locations[5][0] -= width;
354  congruent_locations[5][1] += height;
355 
356  congruent_locations[6][0] -= width;
357  congruent_locations[6][1] -= height;
358 
359  congruent_locations[7][0] += width;
360  congruent_locations[7][1] -= height;
361 
362  // Loop over all other boundary nodes
363  for (unsigned node_b_counter = 0; node_b_counter < boundary_nodes.size(); node_b_counter++)
364  {
365  // Get the index of the current boundary node, and the corresponding node from the mesh
366  unsigned node_b_idx = boundary_nodes[node_b_counter]->GetIndex();
367  Node<2>* p_mesh_node_b = pMesh->GetNode(node_b_idx);
368 
369  if (p_node_a == p_mesh_node_b)
370  {
371  continue;
372  }
373 
374  c_vector<double, 2> node_b_pos = p_mesh_node_b->rGetLocation();
375 
376  // Loop over the eight possible congruent locations and check for coincidence of location
377  for (unsigned pos = 0 ; pos < congruent_locations.size() ; pos++)
378  {
379  if (norm_inf(node_b_pos - congruent_locations[pos]) < mTol)
380  {
381  // Once we find a congruent location, we replace that node in all containing elements
382  std::set<unsigned> containing_elems = p_mesh_node_b->rGetContainingElementIndices();
383 
384  assert(containing_elems.size() > 0);
385 
386  for (std::set<unsigned>::iterator it = containing_elems.begin() ; it != containing_elems.end() ; ++it)
387  {
388  VertexElement<2,2>* p_this_elem = pMesh->GetElement(*it);
389  unsigned local_idx = p_this_elem->GetNodeLocalIndex(p_mesh_node_b->GetIndex());
390 
391  assert(local_idx < UINT_MAX);
392 
393  p_this_elem->AddNode(p_node_a, local_idx);
394  p_this_elem->DeleteNode(local_idx);
395  }
396 
397  // Delete the node_b and remove it from the mesh, and return true to start checking again from scratch
398  pMesh->DeleteNodePriorToReMesh(p_mesh_node_b->GetIndex());
399  pMesh->RemoveDeletedNodes();
400  return true;
401  }
402  }
403  }
404 
405  /*
406  * If we have checked p_node_a against all node_b candidates and have not found a congruent location, p_node_a can
407  * be re-tagged as a non-boundary node, and so will not get checked again next time this function is called.
408  */
409  p_node_a->SetAsBoundaryNode(false);
410 
411  return true;
412 }
413 
414 std::vector<c_vector<double, 2> > VoronoiVertexMeshGenerator::GetInitialPointLocations()
415 {
416  // Create a vector which contains mTotalNumElements spaces
417  std::vector<c_vector<double, 2> > seed_points(mTotalNumElements);
418 
420 
421  // Create the correct number of suitably scaled random numbers
422  for (unsigned point_idx = 0 ; point_idx < mTotalNumElements ; point_idx++)
423  {
424  seed_points[point_idx][0] = p_rand_gen->ranf() * mMultiplierInX;
425  seed_points[point_idx][1] = p_rand_gen->ranf() * mMultiplierInY;
426  }
427 
428  return seed_points;
429 }
430 
431 std::vector<c_vector<double, 2> > VoronoiVertexMeshGenerator::GetElementCentroidsFromMesh()
432 {
433  assert(mpMesh->GetNumElements() == mNumElementsX * mNumElementsY);
434 
435  std::vector<c_vector<double, 2> > element_centroids;
436 
437  // Loop over all elements in the mesh
438  for (unsigned elem_idx = 0; elem_idx < mpMesh->GetNumElements(); elem_idx++)
439  {
440  // Get the current centroid of the element
441  c_vector<double, 2> this_centroid = mpMesh->GetCentroidOfElement(elem_idx);
442 
449  // Account for possible wrap-around in the x-direction
450  if (this_centroid[0] < 0.0)
451  {
452  this_centroid[0] += mMultiplierInX;
453  }
454  else if (this_centroid[0] > mMultiplierInX)
455  {
456  this_centroid[0] -= mMultiplierInX;
457  }
458 
459  // Account for possible wrap-around in the y-direction
460  if (this_centroid[1] < 0.0)
461  {
462  this_centroid[1] += mMultiplierInY;
463  }
464  else if (this_centroid[1] > mMultiplierInY)
465  {
466  this_centroid[1] -= mMultiplierInY;
467  }
468 
469  element_centroids.push_back(this_centroid);
470  }
471 
472  return element_centroids;
473 }
474 
475 void VoronoiVertexMeshGenerator::CreateVoronoiTessellation(std::vector<c_vector<double, 2> >& rSeedLocations)
476 {
477  // Clear the mesh nodes and elements, as they will be replaced in this method
478  std::vector<Node<2>*> nodes;
479  std::vector<VertexElement<2, 2>* > elements;
480 
488  // Datatype is boost::polygon::point_data<int>
489  std::vector<boost_point> points;
490 
491  // Take the locations and map them to integers in a subset of [0, INT_MAX/2] x [0, INT_MAX/2]
492  for (unsigned point_idx = 0 ; point_idx < rSeedLocations.size() ; point_idx++)
493  {
494  // Calculate the correct integer gridpoints
495  int point_x = int( floor( (rSeedLocations[point_idx][0] * mSamplingMultiplier) + 0.5) );
496  int point_y = int( floor( (rSeedLocations[point_idx][1] * mSamplingMultiplier) + 0.5) );
497 
498  points.push_back( boost_point(point_x, point_y) );
499  }
500 
501  // Define offset vectors for the 3 by 3 tessellation
502  std::vector<boost_point> offsets;
503  offsets.push_back( boost_point(-mMaxIntX, mMaxIntY) );
504  offsets.push_back( boost_point( 0, mMaxIntY) );
505  offsets.push_back( boost_point( mMaxIntX, mMaxIntY) );
506  offsets.push_back( boost_point( mMaxIntX, 0) );
507  offsets.push_back( boost_point( mMaxIntX, -mMaxIntY) );
508  offsets.push_back( boost_point( 0, -mMaxIntY) );
509  offsets.push_back( boost_point(-mMaxIntX, -mMaxIntY) );
510  offsets.push_back( boost_point(-mMaxIntX, 0) );
511 
512  // Add all the points in the tessellation to the vector of boost points so in total there are 9 copies of each
513  // seed location, suitably tiled.
514  for (unsigned rep = 0; rep < offsets.size(); rep++)
515  {
516  boost_point offset = offsets[rep];
517  for (unsigned point_idx = 0; point_idx < rSeedLocations.size(); point_idx++)
518  {
519  boost_point new_point = boost_point(points[point_idx].x() + offset.x(), points[point_idx].y() + offset.y());
520  points.push_back(new_point);
521  }
522  }
523 
539  // Construct the Voronoi tessellation of these 9 x mTotalNumElements points
540  voronoi_diagram<double> vd;
541  construct_voronoi(points.begin(), points.end(), &vd);
542 
543  // Loop over the cells in the voronoi diagram to find nodes for our mesh
544  for (voronoi_diagram<double>::const_cell_iterator it = vd.cells().begin();
545  it != vd.cells().end();
546  ++it)
547  {
548  // Get a reference to the current cell
549  const voronoi_diagram<double>::cell_type& cell = *it;
550 
551  // The cells we care about are exactly those whose source_index is less than the size of the locations vector
552  // (i.e. those cells in the central portion of the 3x3 tessellation of source points)
553  if (cell.source_index() < rSeedLocations.size())
554  {
555  // We create a vector of nodes, which will be used to create a MutableElement
556  std::vector<Node<2>*> nodes_this_elem;
557 
558  // Loop over the edges of the current cell
559  const voronoi_diagram<double>::edge_type *edge = cell.incident_edge();
560 
561  do
562  {
563  if (edge->is_primary())
564  {
565  // Calculate the location corresponding to the location of vertex0 of the current edge
566  double x_location = (edge->vertex0()->x()) / mSamplingMultiplier;
567  double y_location = (edge->vertex0()->y()) / mSamplingMultiplier;
568 
569  // Create a node at this location. Default to non-boundary node; this will be updated later
570  Node<2>* p_this_node = new Node<2>(nodes.size(), false, x_location, y_location);
571 
579  unsigned existing_node_idx = UINT_MAX;
580 
581  for (unsigned node_idx = 0; node_idx < nodes.size(); node_idx++)
582  {
583  // Grab the existing node location
584  const c_vector<double, 2>& r_existing_node_location = nodes[node_idx]->rGetLocation();
585 
586  // Equality here is determined entirely on coincidence of position
587  if (fabs(r_existing_node_location[0] - p_this_node->rGetLocation()[0]) < DBL_EPSILON)
588  {
589  if (fabs(r_existing_node_location[1] - p_this_node->rGetLocation()[1]) < DBL_EPSILON)
590  {
591  // If the nodes match, return the existing node index
592  existing_node_idx = node_idx;
593  }
594  }
595  }
596 
597  if (existing_node_idx < UINT_MAX)
598  {
599  // The node was already in nodes vector, and its index is the variable 'existing_node'
600  nodes_this_elem.push_back(nodes[existing_node_idx]);
601  delete p_this_node;
602  }
603  else
604  {
605  // The node does not yet exist - we add it to the nodes vector
606  nodes.push_back(p_this_node);
607  nodes_this_elem.push_back(p_this_node);
608  }
609 
610  // Move to the next edge
611  edge = edge->next();
612  }
613 
614  } while (edge != cell.incident_edge());
615 
616  // Add a new VertexElement to the mElements vector
617  elements.push_back(new VertexElement<2,2>(elements.size(), nodes_this_elem));
618  }
619  }
620 
621  // Loop over the cells in the voronoi diagram to identify boundary nodes
622  for (voronoi_diagram<double>::const_cell_iterator it = vd.cells().begin();
623  it != vd.cells().end();
624  ++it)
625  {
626  // Get a reference to the current cell
627  const voronoi_diagram<double>::cell_type& cell = *it;
628 
629  // The cells we care about are exactly those whose source_index is greater than the size of the locations vector
630  // (i.e. those cells in the outside eight portions of the 3x3 tessellation of source points)
631  if (cell.source_index() >= rSeedLocations.size())
632  {
633  // Loop over the edges of the current cell
634  const voronoi_diagram<double>::edge_type *edge = cell.incident_edge();
635 
636  do
637  {
638  /*
639  * Break out of the do-while if there is an infinite edge; we needn't care about cells at the very edge.
640  */
641  if (edge->is_infinite())
642  {
643  break;
644  }
645 
646  if (edge->is_primary())
647  {
648  c_vector<double, 2> vertex_location;
649 
650  // Get the location of vertex0 of the current edge
651  vertex_location[0] = (edge->vertex0()->x()) / mSamplingMultiplier;
652  vertex_location[1] = (edge->vertex0()->y()) / mSamplingMultiplier;
653 
654  /*
655  * Check whether this location coincides with one of our nodes; if it does, it must be a boundary
656  * node
657  */
658  for (unsigned node_idx = 0 ; node_idx < nodes.size() ; node_idx++)
659  {
660  // Grab the existing node location
661  const c_vector<double, 2>& r_existing_node_location = nodes[node_idx]->rGetLocation();
662 
663  // Equality here is determined entirely on coincidence of position
664  if (fabs(r_existing_node_location[0] - vertex_location[0]) < DBL_EPSILON)
665  {
666  if (fabs(r_existing_node_location[1] - vertex_location[1]) < DBL_EPSILON)
667  {
668  // If the locations match, tag the node as being on the boundary
669  nodes[node_idx]->SetAsBoundaryNode(true);
670  }
671  }
672  }
673  // Move to the next edge
674  edge = edge->next();
675  }
676 
677  } while (edge != cell.incident_edge());
678  }
679  }
680 
681  // Create a new mesh with the current vector of nodes and elements
682  if (mpMesh)
683  {
684  delete mpMesh;
685  }
686  mpMesh = new MutableVertexMesh<2,2>(nodes, elements);
687 }
688 
689 void VoronoiVertexMeshGenerator::ValidateInputAndSetMembers()
690 {
691  // Validate inputs
692  if ((mNumElementsX < 2) || (mNumElementsY < 2))
693  {
694  EXCEPTION("Need at least 2 by 2 cells");
695  }
696 
697  if (mElementTargetArea <= 0.0)
698  {
699  EXCEPTION("Specified target area must be strictly positive");
700  }
701 
702  // The total number of elements requested
703  mTotalNumElements = mNumElementsX * mNumElementsY;
704 
705  // Max of the numbers of elements across and up
706  mMaxNumElems = std::max(mNumElementsX, mNumElementsY);
707 
708  // The multipliers necessary in x and y to ensure all seed points are in [0,1]x[0,1]
709  mMultiplierInX = double(mNumElementsX) / double(mMaxNumElems);
710  mMultiplierInY = double(mNumElementsY) / double(mMaxNumElems);
711 
712  // Floating point representation of INT_MAX/2
713  mSamplingMultiplier = 0.5 * double(INT_MAX);
714 
715  // The max integer that a seed point could be mapped to when discretising
716  mMaxIntX = int( floor( (mMultiplierInX * mSamplingMultiplier) + 0.5 ) );
717  mMaxIntY = int( floor( (mMultiplierInY * mSamplingMultiplier) + 0.5 ) );
718 }
719 
720 void VoronoiVertexMeshGenerator::ValidateSeedLocations(std::vector<c_vector<double, 2> >& rSeedLocations)
721 {
722  unsigned num_seeds = rSeedLocations.size();
723 
724  /*
725  * Seeds at least 1.0 / mSamplingMultiplier apart from each other
726  * are acceptable; the 1.5 term below could just as well be any
727  * number that is strictly greater than 1.0.
728  */
729  double safe_distance = 1.5 / mSamplingMultiplier;
730 
731  // If we find a seed that needs to move position, we will move it and start checking again from the beginning
732  bool recheck = true;
733 
734  while (recheck)
735  {
736  // If no seeds needs moving (as is overwhelmingly likely), we do not need to check again
737  recheck = false;
738 
739  // We check each seed against every other, without double counting
740  for (unsigned seed_idx_one = 0; seed_idx_one < num_seeds; seed_idx_one++)
741  {
742  for (unsigned seed_idx_two = seed_idx_one + 1; seed_idx_two < num_seeds; seed_idx_two++)
743  {
744  // We get the distance between the two seeds currently being considered
745  c_vector<double, 2> one_to_two = rSeedLocations[seed_idx_two] - rSeedLocations[seed_idx_one];
746  double distance = norm_2(one_to_two);
747 
748  if (distance < safe_distance)
749  {
750  // We will now need to re-check
751  recheck = true;
752 
753  /*
754  * If the distance is non-zero, we will move the second seed away along the line joining the two
755  * seeds until it's at the safe distance. If the distance is somehow zero, we will just move the
756  * x-location of the second seed by the safe distance.
757  *
758  * In all cases, we also need to ensure the new location is within the boundary box.
759  */
760  if (distance > DBL_EPSILON)
761  {
762  double multiplier = safe_distance / distance;
763  rSeedLocations[seed_idx_two] += (one_to_two * multiplier);
764 
765  rSeedLocations[seed_idx_two][0] = fmod(rSeedLocations[seed_idx_two][0] + mMultiplierInX, mMultiplierInX);
766  rSeedLocations[seed_idx_two][1] = fmod(rSeedLocations[seed_idx_two][1] + mMultiplierInX, mMultiplierInX);
767  }
768  else
769  {
770  rSeedLocations[seed_idx_two][0] += safe_distance;
771  rSeedLocations[seed_idx_two][0] = fmod(rSeedLocations[seed_idx_two][0], mMultiplierInX);
772  }
773  }
774  }
775 
776  // We will re-check now rather than finishing the outer loop first
777  if (recheck)
778  {
779  break;
780  }
781  }
782  }
783 }
784 
785 std::vector<double> VoronoiVertexMeshGenerator::GetPolygonDistribution()
786 {
787  assert(mpMesh != nullptr);
788 
789  // Number of elements in the mesh
790  unsigned num_elems = mpMesh->GetNumElements();
791 
792  // Store the number of each class of polygons
793  std::vector<unsigned> num_polygons(mMaxExpectedNumSidesPerPolygon - 2, 0);
794 
795  // Container to return the polygon distribution
796  std::vector<double> polygon_dist;
797 
798  // Loop over elements in the mesh to get the number of each class of polygon
799  for (unsigned elem_idx = 0; elem_idx < num_elems; elem_idx++)
800  {
801  unsigned num_nodes_this_elem = mpMesh->GetElement(elem_idx)->GetNumNodes();
802 
803  // All polygons are assumed to have 3, 4, 5, ..., mMaxExpectedNumSidesPerPolygon sides
804  assert(num_nodes_this_elem > 2);
805  assert(num_nodes_this_elem <= mMaxExpectedNumSidesPerPolygon);
806 
807  // Increment correct place in counter - triangles in place 0, squares in 1, etc
808  num_polygons[num_nodes_this_elem - 3]++;
809  }
810 
811  // Loop over the vector of polygon numbers and calculate the distribution vector to return
812  unsigned elems_accounted_for = 0;
813  for (unsigned polygon = 0 ; polygon < num_polygons.size() ; polygon++)
814  {
815  elems_accounted_for += num_polygons[polygon];
816 
817  polygon_dist.push_back(static_cast<double>(num_polygons[polygon]) / static_cast<double>(num_elems));
818 
819  // Only fill the vector of polygon distributions to the point where there are none of higher class in the mesh
820  if (elems_accounted_for == num_elems)
821  {
822  break;
823  }
824  }
825 
826  return polygon_dist;
827 }
828 
829 double VoronoiVertexMeshGenerator::GetAreaCoefficientOfVariation()
830 {
831  assert(mpMesh != nullptr);
832 
833  // Number of elements in the mesh, and check there are at least two
834  unsigned num_elems = mpMesh->GetNumElements();
835  assert(num_elems > 1);
836 
837  double var = 0.0;
838 
839  // Loop over elements in the mesh to get the contributions to the variance
840  for (unsigned elem_idx = 0 ; elem_idx < num_elems ; elem_idx++)
841  {
842  double deviation = mpMesh->GetVolumeOfElement(elem_idx) - mElementTargetArea;
843  var += deviation * deviation;
844  }
845 
846  var /= static_cast<double>(num_elems - 1);
847 
848  return sqrt(var) / mElementTargetArea;
849 }
850 
851 void VoronoiVertexMeshGenerator::RefreshSeedsAndRegenerateMesh()
852 {
853  this->GenerateVoronoiMesh();
854 }
855 
856 void VoronoiVertexMeshGenerator::SetMaxExpectedNumSidesPerPolygon(unsigned maxExpectedNumSidesPerPolygon)
857 {
858  mMaxExpectedNumSidesPerPolygon = maxExpectedNumSidesPerPolygon;
859 }
860 
861 unsigned VoronoiVertexMeshGenerator::GetMaxExpectedNumSidesPerPolygon()
862 {
863  return mMaxExpectedNumSidesPerPolygon;
864 }
865 
866 #endif // BOOST_VERSION >= 105200
867 #if BOOST_VERSION < 105200
868 
874 {
875  EXCEPTION("This is a dummy class. Build with Boost version 1.52 or above for functionality.");
876 }
877 
878 #endif // BOOST_VERSION < 105200
Definition: Node.hpp:58
void SetAsBoundaryNode(bool value=true)
Definition: Node.cpp:127
bool IsBoundaryNode() const
Definition: Node.cpp:164
Node< SPACE_DIM > * GetNode(unsigned index) const
#define EXCEPTION(message)
Definition: Exception.hpp:143
std::set< unsigned > & rGetContainingElementIndices()
Definition: Node.cpp:300
VertexElement< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
Definition: VertexMesh.cpp:628
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
void AddNode(Node< SPACE_DIM > *pNode, const unsigned &rIndex)
unsigned GetNumNodes() const
unsigned GetNumNodes() const
void DeleteNodePriorToReMesh(unsigned index)
static RandomNumberGenerator * Instance()
virtual void ReMesh(VertexElementMap &rElementMap)
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition: Node.cpp:139
unsigned GetIndex() const
unsigned GetNumElements() const
void DeleteNode(const unsigned &rIndex)
unsigned GetIndex() const
Definition: Node.cpp:158
unsigned GetNodeLocalIndex(unsigned globalIndex) const