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