36#include "VoronoiVertexMeshGenerator.hpp"
37#include <boost/make_shared.hpp>
38#include <boost/shared_ptr.hpp>
40#if BOOST_VERSION >= 105200
42using boost::polygon::voronoi_builder;
43using boost::polygon::voronoi_diagram;
44typedef boost::polygon::point_data<int> boost_point;
47 unsigned numElementsY,
48 unsigned numRelaxationSteps,
49 double elementTargetArea)
52 mNumElementsX(numElementsX),
53 mNumElementsY(numElementsY),
54 mNumRelaxationSteps(numRelaxationSteps),
55 mElementTargetArea(elementTargetArea),
57 mMaxExpectedNumSidesPerPolygon(17)
59 this->ValidateInputAndSetMembers();
60 this->GenerateVoronoiMesh();
63void VoronoiVertexMeshGenerator::GenerateVoronoiMesh()
72 std::vector<c_vector<double, 2> > seed_locations = this->GetInitialPointLocations();
73 this->ValidateSeedLocations(seed_locations);
76 this->CreateVoronoiTessellation(seed_locations);
82 for (
unsigned relaxation = 0 ; relaxation < mNumRelaxationSteps ; relaxation++)
84 seed_locations.clear();
85 seed_locations = this->GetElementCentroidsFromMesh();
87 this->CreateVoronoiTessellation(seed_locations);
91 double scale_factor =
double(mMaxNumElems) * sqrt(mElementTargetArea);
92 for (
unsigned node_idx = 0 ; node_idx < mpMesh->GetNumNodes() ; node_idx++)
94 c_vector<double, 2>& node_location = mpMesh->GetNode(node_idx)->rGetModifiableLocation();
95 node_location *= scale_factor;
99boost::shared_ptr<MutableVertexMesh<2,2> > VoronoiVertexMeshGenerator::GetMesh()
104boost::shared_ptr<MutableVertexMesh<2,2> > VoronoiVertexMeshGenerator::GetMeshAfterReMesh()
110boost::shared_ptr<Toroidal2dVertexMesh> VoronoiVertexMeshGenerator::GetToroidalMesh()
127 double width = mNumElementsX * sqrt(mElementTargetArea);
128 double height = mNumElementsY * sqrt(mElementTargetArea);
131 std::vector<Node<2>*> new_nodes(mpMesh->GetNumNodes());
132 std::vector<VertexElement<2,2>*> new_elems(mpMesh->GetNumElements());
135 for (
unsigned node_counter = 0 ; node_counter < mpMesh->GetNumNodes() ; node_counter++)
137 Node<2>* p_node_to_copy = mpMesh->GetNode(node_counter);
140 unsigned copy_index = p_node_to_copy->
GetIndex();
142 const c_vector<double, 2>& copy_location = p_node_to_copy->
rGetLocation();
145 assert(copy_index < mpMesh->GetNumNodes());
148 new_nodes[copy_index] =
new Node<2>(copy_index, copy_location, copy_is_boundary);
152 for (
unsigned elem_counter = 0; elem_counter < mpMesh->GetNumElements(); elem_counter++)
157 unsigned copy_index = p_elem_to_copy->
GetIndex();
158 unsigned copy_num_nodes = p_elem_to_copy->
GetNumNodes();
161 assert(copy_index < mpMesh->GetNumElements());
164 std::vector<Node<2>*> nodes_this_elem;
167 for (
unsigned node_local_idx = 0 ; node_local_idx < copy_num_nodes ; 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]);
194 bool re_check =
true;
198 re_check = this->CheckForCongruentNodes(p_temp_mesh, width, height);
209 for (
unsigned node_counter = 0 ; node_counter < p_temp_mesh->
GetNumNodes() ; node_counter++)
214 unsigned copy_index = p_node_to_copy->
GetIndex();
215 const c_vector<double, 2>& copy_location = p_node_to_copy->
rGetLocation();
221 assert(copy_index < p_temp_mesh->GetNumNodes());
224 new_nodes[copy_index] =
new Node<2>(copy_index, copy_location,
false);
228 for (
unsigned elem_counter = 0; elem_counter < p_temp_mesh->
GetNumElements(); elem_counter++)
233 unsigned copy_index = p_elem_to_copy->
GetIndex();
234 unsigned copy_num_nodes = p_elem_to_copy->
GetNumNodes();
237 assert(copy_index < p_temp_mesh->GetNumElements());
240 std::vector<Node<2>*> nodes_this_elem;
243 for (
unsigned node_local_idx = 0 ; node_local_idx < copy_num_nodes ; node_local_idx++)
247 unsigned local_node_global_idx = p_local_node->
GetIndex();
249 nodes_this_elem.push_back(new_nodes[local_node_global_idx]);
263 mpTorMesh = boost::make_shared<Toroidal2dVertexMesh>(width, height, new_nodes, new_elems);
267 c_vector<double, 2> min_x_y;
268 min_x_y[0] = DBL_MAX;
269 min_x_y[1] = DBL_MAX;
272 for (
unsigned node_idx = 0 ; node_idx < mpTorMesh->GetNumNodes() ; node_idx++)
274 const c_vector<double, 2>& r_this_node_location = mpTorMesh->GetNode(node_idx)->rGetLocation();
276 if (r_this_node_location[0] < min_x_y[0])
278 min_x_y[0] = r_this_node_location[0];
280 if (r_this_node_location[1] < min_x_y[1])
282 min_x_y[1] = r_this_node_location[1];
287 for (
unsigned node_idx = 0 ; node_idx < mpTorMesh->GetNumNodes() ; node_idx++)
289 mpTorMesh->GetNode(node_idx)->rGetModifiableLocation() -= min_x_y;
293 for (
unsigned node_idx = 0 ; node_idx < mpTorMesh->GetNumNodes() ; node_idx++)
295 if (mpTorMesh->GetNode(node_idx)->rGetLocation()[0] >= width)
297 mpTorMesh->GetNode(node_idx)->rGetModifiableLocation()[0] -= width;
299 if (mpTorMesh->GetNode(node_idx)->rGetLocation()[1] >= height)
301 mpTorMesh->GetNode(node_idx)->rGetModifiableLocation()[1] -= height;
308bool VoronoiVertexMeshGenerator::CheckForCongruentNodes(
MutableVertexMesh<2,2>* pMesh,
double width,
double height)
311 std::vector<Node<2>*> boundary_nodes;
312 for (
unsigned node_idx = 0 ; node_idx < pMesh->
GetNumNodes() ; node_idx++)
318 boundary_nodes.push_back(p_node);
323 if (boundary_nodes.empty())
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);
333 congruent_locations[0][0] += width;
335 congruent_locations[1][1] += height;
337 congruent_locations[2][0] -= width;
339 congruent_locations[3][1] -= height;
341 congruent_locations[4][0] += width;
342 congruent_locations[4][1] += height;
344 congruent_locations[5][0] -= width;
345 congruent_locations[5][1] += height;
347 congruent_locations[6][0] -= width;
348 congruent_locations[6][1] -= height;
350 congruent_locations[7][0] += width;
351 congruent_locations[7][1] -= height;
354 for (
unsigned node_b_counter = 0; node_b_counter < boundary_nodes.size(); node_b_counter++)
357 unsigned node_b_idx = boundary_nodes[node_b_counter]->GetIndex();
360 if (p_node_a == p_mesh_node_b)
365 c_vector<double, 2> node_b_pos = p_mesh_node_b->
rGetLocation();
368 for (
unsigned pos = 0 ; pos < congruent_locations.size() ; pos++)
370 if (norm_inf(node_b_pos - congruent_locations[pos]) < mTol)
375 assert(containing_elems.size() > 0);
377 for (std::set<unsigned>::iterator it = containing_elems.begin() ; it != containing_elems.end() ; ++it)
382 assert(local_idx < UINT_MAX);
384 p_this_elem->
AddNode(p_node_a, local_idx);
405std::vector<c_vector<double, 2> > VoronoiVertexMeshGenerator::GetInitialPointLocations()
408 std::vector<c_vector<double, 2> > seed_points(mTotalNumElements);
413 for (
unsigned point_idx = 0 ; point_idx < mTotalNumElements ; point_idx++)
415 seed_points[point_idx][0] = p_rand_gen->
ranf() * mMultiplierInX;
416 seed_points[point_idx][1] = p_rand_gen->
ranf() * mMultiplierInY;
422std::vector<c_vector<double, 2> > VoronoiVertexMeshGenerator::GetElementCentroidsFromMesh()
424 assert(mpMesh->GetNumElements() == mNumElementsX * mNumElementsY);
426 std::vector<c_vector<double, 2> > element_centroids;
429 for (
unsigned elem_idx = 0; elem_idx < mpMesh->GetNumElements(); elem_idx++)
432 c_vector<double, 2> this_centroid = mpMesh->GetCentroidOfElement(elem_idx);
441 if (this_centroid[0] < 0.0)
443 this_centroid[0] += mMultiplierInX;
445 else if (this_centroid[0] > mMultiplierInX)
447 this_centroid[0] -= mMultiplierInX;
451 if (this_centroid[1] < 0.0)
453 this_centroid[1] += mMultiplierInY;
455 else if (this_centroid[1] > mMultiplierInY)
457 this_centroid[1] -= mMultiplierInY;
460 element_centroids.push_back(this_centroid);
463 return element_centroids;
466void VoronoiVertexMeshGenerator::CreateVoronoiTessellation(std::vector<c_vector<double, 2> >& rSeedLocations)
469 std::vector<Node<2>*> nodes;
470 std::vector<VertexElement<2, 2>* > elements;
480 std::vector<boost_point> points;
483 for (
unsigned point_idx = 0 ; point_idx < rSeedLocations.size() ; point_idx++)
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) );
489 points.push_back( boost_point(point_x, point_y) );
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) );
505 for (
unsigned rep = 0; rep < offsets.size(); rep++)
507 boost_point offset = offsets[rep];
508 for (
unsigned point_idx = 0; point_idx < rSeedLocations.size(); point_idx++)
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);
531 voronoi_diagram<double> vd;
532 construct_voronoi(points.begin(), points.end(), &vd);
535 for (voronoi_diagram<double>::const_cell_iterator it = vd.cells().begin();
536 it != vd.cells().end();
540 const voronoi_diagram<double>::cell_type& cell = *it;
544 if (cell.source_index() < rSeedLocations.size())
547 std::vector<Node<2>*> nodes_this_elem;
550 const voronoi_diagram<double>::edge_type *edge = cell.incident_edge();
554 if (edge->is_primary())
557 double x_location = (edge->vertex0()->x()) / mSamplingMultiplier;
558 double y_location = (edge->vertex0()->y()) / mSamplingMultiplier;
561 Node<2>* p_this_node =
new Node<2>(nodes.size(),
false, x_location, y_location);
570 unsigned existing_node_idx = UINT_MAX;
572 for (
unsigned node_idx = 0; node_idx < nodes.size(); node_idx++)
575 const c_vector<double, 2>& r_existing_node_location = nodes[node_idx]->rGetLocation();
578 if (fabs(r_existing_node_location[0] - p_this_node->
rGetLocation()[0]) < DBL_EPSILON)
580 if (fabs(r_existing_node_location[1] - p_this_node->
rGetLocation()[1]) < DBL_EPSILON)
583 existing_node_idx = node_idx;
588 if (existing_node_idx < UINT_MAX)
591 nodes_this_elem.push_back(nodes[existing_node_idx]);
597 nodes.push_back(p_this_node);
598 nodes_this_elem.push_back(p_this_node);
605 }
while (edge != cell.incident_edge());
613 for (voronoi_diagram<double>::const_cell_iterator it = vd.cells().begin();
614 it != vd.cells().end();
618 const voronoi_diagram<double>::cell_type& cell = *it;
622 if (cell.source_index() >= rSeedLocations.size())
625 const voronoi_diagram<double>::edge_type *edge = cell.incident_edge();
632 if (edge->is_infinite())
637 if (edge->is_primary())
639 c_vector<double, 2> vertex_location;
642 vertex_location[0] = (edge->vertex0()->x()) / mSamplingMultiplier;
643 vertex_location[1] = (edge->vertex0()->y()) / mSamplingMultiplier;
649 for (
unsigned node_idx = 0 ; node_idx < nodes.size() ; node_idx++)
652 const c_vector<double, 2>& r_existing_node_location = nodes[node_idx]->
rGetLocation();
655 if (fabs(r_existing_node_location[0] - vertex_location[0]) < DBL_EPSILON)
657 if (fabs(r_existing_node_location[1] - vertex_location[1]) < DBL_EPSILON)
660 nodes[node_idx]->SetAsBoundaryNode(
true);
668 }
while (edge != cell.incident_edge());
677 mpMesh = boost::make_shared<MutableVertexMesh<2,2> >(nodes, elements);
680void VoronoiVertexMeshGenerator::ValidateInputAndSetMembers()
683 if ((mNumElementsX < 2) || (mNumElementsY < 2))
688 if (mElementTargetArea <= 0.0)
690 EXCEPTION(
"Specified target area must be strictly positive");
694 mTotalNumElements = mNumElementsX * mNumElementsY;
697 mMaxNumElems = std::max(mNumElementsX, mNumElementsY);
700 mMultiplierInX =
double(mNumElementsX) /
double(mMaxNumElems);
701 mMultiplierInY =
double(mNumElementsY) /
double(mMaxNumElems);
704 mSamplingMultiplier = 0.5 *
double(INT_MAX);
707 mMaxIntX = int( floor( (mMultiplierInX * mSamplingMultiplier) + 0.5 ) );
708 mMaxIntY = int( floor( (mMultiplierInY * mSamplingMultiplier) + 0.5 ) );
711void VoronoiVertexMeshGenerator::ValidateSeedLocations(std::vector<c_vector<double, 2> >& rSeedLocations)
713 unsigned num_seeds = rSeedLocations.size();
720 double safe_distance = 1.5 / mSamplingMultiplier;
731 for (
unsigned seed_idx_one = 0; seed_idx_one < num_seeds; seed_idx_one++)
733 for (
unsigned seed_idx_two = seed_idx_one + 1; seed_idx_two < num_seeds; seed_idx_two++)
736 c_vector<double, 2> one_to_two = rSeedLocations[seed_idx_two] - rSeedLocations[seed_idx_one];
737 double distance = norm_2(one_to_two);
739 if (distance < safe_distance)
751 if (distance > DBL_EPSILON)
753 double multiplier = safe_distance / distance;
754 rSeedLocations[seed_idx_two] += (one_to_two * multiplier);
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);
761 rSeedLocations[seed_idx_two][0] += safe_distance;
762 rSeedLocations[seed_idx_two][0] = fmod(rSeedLocations[seed_idx_two][0], mMultiplierInX);
776std::vector<double> VoronoiVertexMeshGenerator::GetPolygonDistribution()
778 assert(mpMesh !=
nullptr);
781 unsigned num_elems = mpMesh->GetNumElements();
784 std::vector<unsigned> num_polygons(mMaxExpectedNumSidesPerPolygon - 2, 0);
787 std::vector<double> polygon_dist;
790 for (
unsigned elem_idx = 0; elem_idx < num_elems; elem_idx++)
792 unsigned num_nodes_this_elem = mpMesh->GetElement(elem_idx)->GetNumNodes();
795 assert(num_nodes_this_elem > 2);
796 assert(num_nodes_this_elem <= mMaxExpectedNumSidesPerPolygon);
799 num_polygons[num_nodes_this_elem - 3]++;
803 unsigned elems_accounted_for = 0;
804 for (
unsigned polygon = 0 ; polygon < num_polygons.size() ; polygon++)
806 elems_accounted_for += num_polygons[polygon];
808 polygon_dist.push_back(
static_cast<double>(num_polygons[polygon]) /
static_cast<double>(num_elems));
811 if (elems_accounted_for == num_elems)
820double VoronoiVertexMeshGenerator::GetAreaCoefficientOfVariation()
822 assert(mpMesh !=
nullptr);
825 unsigned num_elems = mpMesh->GetNumElements();
826 assert(num_elems > 1);
831 for (
unsigned elem_idx = 0 ; elem_idx < num_elems ; elem_idx++)
833 double deviation = mpMesh->GetVolumeOfElement(elem_idx) - mElementTargetArea;
834 var += deviation * deviation;
837 var /=
static_cast<double>(num_elems - 1);
839 return sqrt(var) / mElementTargetArea;
842void VoronoiVertexMeshGenerator::RefreshSeedsAndRegenerateMesh()
844 this->GenerateVoronoiMesh();
847void VoronoiVertexMeshGenerator::SetMaxExpectedNumSidesPerPolygon(
unsigned maxExpectedNumSidesPerPolygon)
849 mMaxExpectedNumSidesPerPolygon = maxExpectedNumSidesPerPolygon;
852unsigned VoronoiVertexMeshGenerator::GetMaxExpectedNumSidesPerPolygon()
854 return mMaxExpectedNumSidesPerPolygon;
858#if BOOST_VERSION < 105200
866 EXCEPTION(
"This is a dummy class. Build with Boost version 1.52 or above for functionality.");
#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)
void RemoveDeletedNodes()
unsigned GetNumNodes() const
void DeleteNodePriorToReMesh(unsigned index)
unsigned GetNumElements() const
std::set< unsigned > & rGetContainingElementIndices()
const c_vector< double, SPACE_DIM > & rGetLocation() const
bool IsBoundaryNode() const
unsigned GetIndex() const
void SetAsBoundaryNode(bool value=true)
static RandomNumberGenerator * Instance()
VertexElement< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
VoronoiVertexMeshGenerator()