Chaste  Build::
Part.cpp
1 /*
2 
3 Copyright (c) 2005-2016, 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 <math.h>
37 #define _BACKWARD_BACKWARD_WARNING_H 1 //Cut out the vtk deprecated warning
38 #include "vtkCleanPolyData.h"
39 #include "vtkSelectEnclosedPoints.h"
40 #include <vtkVersion.h>
41 #include <boost/random.hpp>
42 #include <boost/generator_iterator.hpp>
43 #include "VesselNode.hpp"
44 #include "VesselSegment.hpp"
45 #include "Vessel.hpp"
46 #include "VesselSurfaceGenerator.hpp"
47 #include "Part.hpp"
48 #include "GeometryTools.hpp"
49 #include "BaseUnits.hpp"
50 
51 template<unsigned DIM>
53  mFacets(),
54  mVtkPart(vtkSmartPointer<vtkPolyData>()),
55  mHoleMarkers(),
56  mRegionMarkers(),
57  mReferenceLength(BaseUnits::Instance()->GetReferenceLengthScale()),
58  mVtkIsUpToDate(false)
59 {
60 }
61 
62 template<unsigned DIM>
63 boost::shared_ptr<Part<DIM> > Part<DIM>::Create()
64 {
65  MAKE_PTR(Part<DIM>, pSelf);
66  return pSelf;
67 }
68 
69 template<unsigned DIM>
71 {
72 
73 }
74 
75 template<unsigned DIM>
76 boost::shared_ptr<Polygon<DIM> > Part<DIM>::AddCircle(units::quantity<unit::length> radius,
77  DimensionalChastePoint<DIM> centre, unsigned numSegments)
78 {
79  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > vertices;
80  double seg_angle = 2.0 * M_PI / double(numSegments);
81  for (unsigned idx = 0; idx < numSegments; idx++)
82  {
83  double angle = seg_angle * double(idx);
84  double x = (radius * std::cos(angle) + centre.GetLocation(mReferenceLength)[0]*centre.GetReferenceLengthScale())/mReferenceLength;
85  double y = (radius * std::sin(angle) + centre.GetLocation(mReferenceLength)[1]*centre.GetReferenceLengthScale())/mReferenceLength;
86  if(DIM==3)
87  {
89  }
90  else
91  {
92  vertices.push_back(DimensionalChastePoint<DIM>::Create(x, y, 0.0, mReferenceLength));
93  }
94  }
95  mVtkIsUpToDate = false;
96  return AddPolygon(vertices);
97 }
98 
99 template<unsigned DIM>
100 void Part<DIM>::AddCylinder(units::quantity<unit::length> radius,
101  units::quantity<unit::length> depth,
103  unsigned numSegments)
104 {
105  boost::shared_ptr<Polygon<DIM> > p_circle = AddCircle(radius, centre, numSegments);
106  Extrude(p_circle, depth);
107  mVtkIsUpToDate = false;
108 }
109 
110 template<unsigned DIM>
111 void Part<DIM>::AddCuboid(units::quantity<unit::length> sizeX,
112  units::quantity<unit::length> sizeY,
113  units::quantity<unit::length> sizeZ,
115 {
116  boost::shared_ptr<Polygon<DIM> > p_rectangle = AddRectangle(sizeX, sizeY, origin);
117  Extrude(p_rectangle, sizeZ);
118  mVtkIsUpToDate = false;
119 }
120 
121 template<unsigned DIM>
123 {
124  mHoleMarkers.push_back(hole);
125 }
126 
127 template<unsigned DIM>
128 boost::shared_ptr<Polygon<DIM> > Part<DIM>::AddPolygon(std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > vertices, bool newFacet,
129  boost::shared_ptr<Facet<DIM> > pFacet)
130 {
131  boost::shared_ptr<Polygon<DIM> > p_polygon = Polygon<DIM>::Create(vertices);
132  AddPolygon(p_polygon, newFacet, pFacet);
133  mVtkIsUpToDate = false;
134  return p_polygon;
135 }
136 
137 template<unsigned DIM>
138 boost::shared_ptr<Polygon<DIM> > Part<DIM>::AddPolygon(boost::shared_ptr<Polygon<DIM> > pPolygon, bool newFacet,
139  boost::shared_ptr<Facet<DIM> > pFacet)
140 {
141  if (!pFacet)
142  {
143  if (mFacets.size() == 0 || newFacet)
144  {
145  mFacets.push_back(Facet<DIM>::Create(pPolygon));
146  }
147  else
148  {
149  mFacets[0]->AddPolygon(pPolygon);
150  }
151  }
152  else
153  {
154  pFacet->AddPolygon(pPolygon);
155  }
156  mVtkIsUpToDate = false;
157  return pPolygon;
158 }
159 
160 template<unsigned DIM>
161 boost::shared_ptr<Polygon<DIM> > Part<DIM>::AddRectangle(units::quantity<unit::length> sizeX,
162  units::quantity<unit::length> sizeY,
164 {
165  c_vector<double, DIM> dimensionless_origin = origin.GetLocation(mReferenceLength);
166  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > vertices;
167  if(DIM==3)
168  {
169  vertices.push_back(DimensionalChastePoint<DIM>::Create(dimensionless_origin[0], dimensionless_origin[1], dimensionless_origin[2], mReferenceLength));
170  vertices.push_back(DimensionalChastePoint<DIM>::Create(dimensionless_origin[0] + sizeX/mReferenceLength,
171  dimensionless_origin[1], dimensionless_origin[2],mReferenceLength));
172  vertices.push_back(DimensionalChastePoint<DIM>::Create(dimensionless_origin[0] + sizeX/mReferenceLength,
173  dimensionless_origin[1] + sizeY/mReferenceLength, dimensionless_origin[2], mReferenceLength));
174  vertices.push_back(DimensionalChastePoint<DIM>::Create(dimensionless_origin[0], dimensionless_origin[1] + sizeY/mReferenceLength, dimensionless_origin[2],mReferenceLength));
175  }
176  else
177  {
178  vertices.push_back(DimensionalChastePoint<DIM>::Create(dimensionless_origin[0], dimensionless_origin[1], 0.0, mReferenceLength));
179  vertices.push_back(DimensionalChastePoint<DIM>::Create(dimensionless_origin[0] + sizeX/mReferenceLength, dimensionless_origin[1], 0.0, mReferenceLength));
180  vertices.push_back(DimensionalChastePoint<DIM>::Create(dimensionless_origin[0] + sizeX/mReferenceLength,
181  dimensionless_origin[1] + sizeY/mReferenceLength, 0.0, mReferenceLength));
182  vertices.push_back(DimensionalChastePoint<DIM>::Create(dimensionless_origin[0], dimensionless_origin[1] + sizeY/mReferenceLength, 0.0, mReferenceLength));
183  }
184  mVtkIsUpToDate = false;
185  return AddPolygon(vertices);
186 }
187 
188 template<unsigned DIM>
189 void Part<DIM>::AddVesselNetwork(boost::shared_ptr<VesselNetwork<DIM> > pVesselNetwork, bool surface)
190 {
191  if (!surface)
192  {
193  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > vertices;
194  std::vector<boost::shared_ptr<VesselNode<DIM> > > nodes = pVesselNetwork->GetNodes();
195  for (unsigned idx = 0; idx < nodes.size(); idx++)
196  {
197  units::quantity<unit::length> length_Scale = nodes[idx]->rGetLocation().GetReferenceLengthScale();
198  vertices.push_back(DimensionalChastePoint<DIM>::Create(nodes[idx]->rGetLocation().GetLocation(length_Scale), length_Scale));
199  }
200 
201  // If vertices lie on any existing facets add the vertex to the facet
202  for (unsigned kdx = 0; kdx < vertices.size(); kdx++)
203  {
204  for (unsigned idx = 0; idx < mFacets.size(); idx++)
205  {
206  if(mFacets[idx]->ContainsPoint(*(vertices[kdx])))
207  {
208  mFacets[idx]->AddPolygon(Polygon<DIM>::Create(vertices[kdx]));
209  }
210  }
211  }
212 
213  // Create polygons and facets for each vessel
214  std::vector<boost::shared_ptr<Vessel<DIM> > > vessels = pVesselNetwork->GetVessels();
215  for (unsigned idx = 0; idx < vessels.size(); idx++)
216  {
217  std::vector<boost::shared_ptr<VesselSegment<DIM> > > segments = vessels[idx]->GetSegments();
218  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > segment_vertices;
219  for (unsigned jdx = 0; jdx < segments.size(); jdx++)
220  {
221  unsigned node0_index = pVesselNetwork->GetNodeIndex(segments[jdx]->GetNode(0));
222  unsigned node1_index = pVesselNetwork->GetNodeIndex(segments[jdx]->GetNode(1));
223  if (jdx == 0)
224  {
225  segment_vertices.push_back(vertices[node0_index]);
226  }
227  segment_vertices.push_back(vertices[node1_index]);
228  }
229  boost::shared_ptr<Polygon<DIM> > p_polygon = Polygon<DIM>::Create(segment_vertices);
230  mFacets.push_back(Facet<DIM>::Create(p_polygon));
231  }
232  }
233  else
234  {
235  // Add any polygons on existing facets to the facet
236  VesselSurfaceGenerator<DIM> generator(pVesselNetwork);
237  std::vector<boost::shared_ptr<Polygon<DIM> > > polygons = generator.GetSurfacePolygons();
238  std::vector<bool> polygon_on_facet;
239 
240  for (unsigned idx = 0; idx < polygons.size(); idx++)
241  {
242  bool on_facet = false;
243  DimensionalChastePoint<DIM> poly_centroid = polygons[idx]->GetCentroid();
244  for (unsigned jdx = 0; jdx < mFacets.size(); jdx++)
245  {
246  if (mFacets[jdx]->ContainsPoint(poly_centroid))
247  {
248  on_facet = true;
249  mFacets[jdx]->AddPolygon(polygons[idx]);
250  }
251  }
252  polygon_on_facet.push_back(on_facet);
253  }
254 
255  // Create polygons and facets for each vessel
256  for (unsigned idx = 0; idx < polygons.size(); idx++)
257  {
258  if (!polygon_on_facet[idx])
259  {
260  mFacets.push_back(Facet<DIM>::Create(polygons[idx]));
261  }
262  }
263 
264  std::vector<DimensionalChastePoint<DIM> > hole_locations = generator.GetHoles();
265  for(unsigned idx=0; idx<hole_locations.size(); idx++)
266  {
267  AddHoleMarker(hole_locations[idx]);
268  }
269  }
270  mVtkIsUpToDate = false;
271 }
272 
273 template<unsigned DIM>
274 void Part<DIM>::Extrude(boost::shared_ptr<Polygon<DIM> > pPolygon, units::quantity<unit::length> depth)
275 {
276  if(DIM==2)
277  {
278  EXCEPTION("Only parts in 3D space can be extruded.");
279  }
280  // Loop through the vertices and create new ones at the offset depth
281  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > original_vertices = pPolygon->GetVertices();
282  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > new_vertices;
283  for (unsigned idx = 0; idx < original_vertices.size(); idx++)
284  {
285  c_vector<double, DIM> vertex_location = original_vertices[idx]->GetLocation(mReferenceLength);
286  if(DIM==2)
287  {
288  new_vertices.push_back(DimensionalChastePoint<DIM>::Create(vertex_location[0], vertex_location[1], depth/mReferenceLength, mReferenceLength));
289  }
290  else
291  {
292  new_vertices.push_back(DimensionalChastePoint<DIM>::Create(vertex_location[0], vertex_location[1],
293  vertex_location[2] + depth/mReferenceLength, mReferenceLength));
294  }
295  }
296 
297  // Every straight edge is now a planar face, with 3 new edges ordered in CCW
298  for (unsigned idx = 0; idx < original_vertices.size(); idx++)
299  {
300  unsigned index2;
301  if (idx != original_vertices.size() - 1)
302  {
303  index2 = idx + 1;
304  }
305  else
306  {
307  index2 = 0;
308  }
309  boost::shared_ptr<Polygon<DIM> > p_polygon = Polygon<DIM>::Create(original_vertices[idx]);
310  p_polygon->AddVertex(original_vertices[index2]);
311  p_polygon->AddVertex(new_vertices[index2]);
312  p_polygon->AddVertex(new_vertices[idx]);
313  mFacets.push_back(Facet<DIM>::Create(p_polygon));
314  }
315 
316  // Close the lid
317  boost::shared_ptr<Polygon<DIM> > p_polygon = Polygon<DIM>::Create(new_vertices);
318  mFacets.push_back(Facet<DIM>::Create(p_polygon));
319  mVtkIsUpToDate = false;
320 }
321 
322 template <unsigned DIM>
323 void Part<DIM>::BooleanWithNetwork(boost::shared_ptr<VesselNetwork<DIM> > pVesselNetwork)
324 {
325  // Remove any vessel with both nodes outside the domain
326  std::vector<boost::shared_ptr<Vessel<DIM> > > vessels = pVesselNetwork->GetVessels();
327  for(unsigned idx=0;idx<vessels.size();idx++)
328  {
329  if(!IsPointInPart(vessels[idx]->GetStartNode()->rGetLocation()) &&
330  !IsPointInPart(vessels[idx]->GetEndNode()->rGetLocation()))
331  {
332  pVesselNetwork->RemoveVessel(vessels[idx], true);
333  }
334  }
335  pVesselNetwork->UpdateAll();
336 }
337 
338 template<unsigned DIM>
339 std::vector<DimensionalChastePoint<DIM> > Part<DIM>::GetHoleMarkers()
340 {
341  return mHoleMarkers;
342 }
343 
344 template<unsigned DIM>
345 units::quantity<unit::length> Part<DIM>::GetReferenceLengthScale()
346 {
347  return mReferenceLength;
348 }
349 
350 template<unsigned DIM>
351 boost::shared_ptr<Facet<DIM> > Part<DIM>::GetFacet(const DimensionalChastePoint<DIM>& location)
352 {
353  for(unsigned idx=0; idx<mFacets.size(); idx++)
354  {
355  if(mFacets[idx]->ContainsPoint(location))
356  {
357  return mFacets[idx];
358  }
359  }
360  EXCEPTION("No facet found at input location");
361 }
362 
363 template<unsigned DIM>
364 std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > Part<DIM>::GetVertices()
365 {
366  std::vector<boost::shared_ptr<Polygon<DIM> > > polygons = GetPolygons();
367  std::set<boost::shared_ptr<DimensionalChastePoint<DIM> > > unique_vertices;
368  for (unsigned idx = 0; idx < polygons.size(); idx++)
369  {
370  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > polygon_vertices = polygons[idx]->GetVertices();
371  std::copy(polygon_vertices.begin(), polygon_vertices.end(),
372  std::inserter(unique_vertices, unique_vertices.end()));
373  }
374 
375  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > vertices;
376  vertices.insert(vertices.end(), unique_vertices.begin(), unique_vertices.end());
377 
378  for (unsigned idx = 0; idx < vertices.size(); idx++)
379  {
380  vertices[idx]->SetIndex(idx);
381  }
382  return vertices;
383 }
384 
385 template<unsigned DIM>
386 std::vector<unsigned> Part<DIM>::GetContainingGridIndices(unsigned num_x, unsigned num_y, unsigned num_z, units::quantity<unit::length> spacing)
387 {
388  double scaled_spacing = spacing/mReferenceLength;
389  std::vector<unsigned> location_indices;
390  for(unsigned kdx=0; kdx<num_z; kdx++)
391  {
392  for(unsigned jdx=0; jdx<num_y; jdx++)
393  {
394  for(unsigned idx=0; idx<num_x; idx++)
395  {
396  DimensionalChastePoint<DIM> location(double(idx) * scaled_spacing, double(jdx) * scaled_spacing, double(kdx) * scaled_spacing, mReferenceLength);
397  unsigned index = idx + num_x * jdx + num_x * num_y * kdx;
398  if(IsPointInPart(location))
399  {
400  location_indices.push_back(index);
401  }
402  }
403  }
404  }
405  return location_indices;
406 }
407 
408 template<unsigned DIM>
409 std::vector<DimensionalChastePoint<DIM> > Part<DIM>::GetVertexLocations()
410 {
411  std::vector<DimensionalChastePoint<DIM> > vertex_locs;
412  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > vertices = GetVertices();
413  for (unsigned idx = 0; idx < vertices.size(); idx++)
414  {
415  vertex_locs.push_back(*vertices[idx]);
416  }
417  return vertex_locs;
418 }
419 
420 template<unsigned DIM>
421 std::vector<boost::shared_ptr<Polygon<DIM> > > Part<DIM>::GetPolygons()
422 {
423  std::vector<boost::shared_ptr<Polygon<DIM> > > polygons;
424  for (unsigned idx = 0; idx < mFacets.size(); idx++)
425  {
426  std::vector<boost::shared_ptr<Polygon<DIM> > > facet_polygons = mFacets[idx]->GetPolygons();
427  polygons.insert(polygons.end(), facet_polygons.begin(), facet_polygons.end());
428  }
429  return polygons;
430 }
431 
432 template<unsigned DIM>
433 std::vector<units::quantity<unit::length> > Part<DIM>::GetBoundingBox()
434 {
435  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > vertices = GetVertices();
436  c_vector<double, 6> box;
437 
438  for (unsigned idx = 0; idx < vertices.size(); idx++)
439  {
440  c_vector<double, DIM> vertex_location = vertices[idx]->GetLocation(mReferenceLength);
441  for (unsigned jdx = 0; jdx < DIM; jdx++)
442  {
443  if (idx == 0)
444  {
445  box[2 * jdx] = vertex_location[jdx];
446  box[2 * jdx + 1] = vertex_location[jdx];
447  }
448  else
449  {
450  if (vertex_location[jdx] < box[2 * jdx])
451  {
452  box[2 * jdx] = vertex_location[jdx];
453  }
454  if (vertex_location[jdx] > box[2 * jdx + 1])
455  {
456  box[2 * jdx + 1] = vertex_location[jdx];
457  }
458  }
459  }
460  }
461 
462  std::vector<units::quantity<unit::length> > box_vector(6, 0.0*unit::metres);
463  for(unsigned idx=0; idx<6; idx++)
464  {
465  box_vector[idx] = box[idx] * mReferenceLength;
466  }
467 
468  return box_vector;
469 }
470 
471 template<unsigned DIM>
472 std::vector<boost::shared_ptr<Facet<DIM> > > Part<DIM>::GetFacets()
473 {
474  return mFacets;
475 }
476 
477 template<unsigned DIM>
478 std::vector<std::pair<unsigned, unsigned> > Part<DIM>::GetSegmentIndices()
479 {
480  // Make sure the vertex indexes are up-to-date.
481  GetVertices();
482 
483  std::vector<std::pair<unsigned, unsigned> > indexes;
484  std::vector<boost::shared_ptr<Polygon<DIM> > > polygons = mFacets[0]->GetPolygons();
485 
486  for (unsigned idx = 0; idx < polygons.size(); idx++)
487  {
488  if (polygons[idx]->GetVertices().size() == 2)
489  {
490  indexes.push_back(
491  std::pair<unsigned, unsigned>(polygons[idx]->GetVertices()[0]->GetIndex(),
492  polygons[idx]->GetVertices()[1]->GetIndex()));
493  }
494  else if (polygons[idx]->GetVertices().size() > 1)
495  {
496  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > vertices = polygons[idx]->GetVertices();
497  for (unsigned jdx = 0; jdx < vertices.size() - 1; jdx++)
498  {
499  indexes.push_back(
500  std::pair<unsigned, unsigned>(vertices[jdx]->GetIndex(), vertices[jdx + 1]->GetIndex()));
501  }
502  indexes.push_back(
503  std::pair<unsigned, unsigned>(vertices[vertices.size() - 1]->GetIndex(), vertices[0]->GetIndex()));
504 
505  }
506  }
507  return indexes;
508 }
509 
510 template<unsigned DIM>
511 vtkSmartPointer<vtkPolyData> Part<DIM>::GetVtk()
512 {
513  if(mVtkIsUpToDate)
514  {
515  return mVtkPart;
516  }
517 
518  vtkSmartPointer<vtkPolyData> p_part_data = vtkSmartPointer<vtkPolyData>::New();
519  vtkSmartPointer<vtkPoints> p_vertices = vtkSmartPointer<vtkPoints>::New();
520 
521  p_part_data->Allocate(1, 1);
522  // Loop through each polygon, collect the vertices and set correct point ids
523 
524  std::vector<boost::shared_ptr<Polygon<DIM> > > polygons = GetPolygons();
525  unsigned vert_counter = 0;
526  for (vtkIdType idx = 0; idx < vtkIdType(polygons.size()); idx++)
527  {
528  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > vertices = polygons[idx]->GetVertices();
529  vtkSmartPointer<vtkPolygon> p_polygon = vtkSmartPointer<vtkPolygon>::New();
530  p_polygon->GetPointIds()->SetNumberOfIds(vertices.size());
531  for (vtkIdType jdx = 0; jdx < vtkIdType(vertices.size()); jdx++)
532  {
533  c_vector<double, DIM> vertex_location = vertices[jdx]->GetLocation(mReferenceLength);
534  if(DIM==3)
535  {
536  p_vertices->InsertNextPoint(vertex_location[0], vertex_location[1], vertex_location[2]);
537  }
538  else
539  {
540  p_vertices->InsertNextPoint(vertex_location[0], vertex_location[1], 0.0);
541  }
542 
543  p_polygon->GetPointIds()->SetId(jdx, vert_counter);
544  vert_counter++;
545  }
546  p_part_data->InsertNextCell(p_polygon->GetCellType(), p_polygon->GetPointIds());
547  }
548  p_part_data->SetPoints(p_vertices);
549 
550  vtkSmartPointer<vtkCleanPolyData> p_clean_data = vtkSmartPointer<vtkCleanPolyData>::New();
551  #if VTK_MAJOR_VERSION <= 5
552  p_clean_data->SetInput(p_part_data);
553  #else
554  p_clean_data->SetInputData(p_part_data);
555  #endif
556  p_clean_data->Update();
557 
558  mVtkPart = p_clean_data->GetOutput();
559  mVtkIsUpToDate = true;
560  return mVtkPart;
561 }
562 
563 template<unsigned DIM>
565 {
566  vtkSmartPointer<vtkPolyData> p_part = GetVtk();
567  vtkSmartPointer<vtkPoints> p_points = vtkSmartPointer<vtkPoints>::New();
568 
569  if(DIM==3)
570  {
571  p_points->InsertNextPoint(location.GetLocation(mReferenceLength)[0], location.GetLocation(mReferenceLength)[1], location.GetLocation(mReferenceLength)[2]);
572  }
573  else
574  {
575  p_points->InsertNextPoint(location.GetLocation(mReferenceLength)[0], location.GetLocation(mReferenceLength)[1], 0.0);
576  }
577 
578  vtkSmartPointer<vtkPolyData> p_point_data = vtkSmartPointer<vtkPolyData>::New();
579  p_point_data->SetPoints(p_points);
580 
581  //Points inside test
582  vtkSmartPointer<vtkSelectEnclosedPoints> selectEnclosedPoints = vtkSmartPointer<vtkSelectEnclosedPoints>::New();
583  #if VTK_MAJOR_VERSION <= 5
584  selectEnclosedPoints->SetInput(p_point_data);
585  #else
586  selectEnclosedPoints->SetInputData(p_point_data);
587  #endif
588  #if VTK_MAJOR_VERSION <= 5
589  selectEnclosedPoints->SetSurface(p_part);
590  #else
591  selectEnclosedPoints->SetSurfaceData(p_part);
592  #endif
593  selectEnclosedPoints->Update();
594 
595  return selectEnclosedPoints->IsInside(0);
596 }
597 
598 template<unsigned DIM>
599 std::vector<bool> Part<DIM>::IsPointInPart(const std::vector<DimensionalChastePoint<DIM> >& rLocations)
600 {
601  vtkSmartPointer<vtkPolyData> p_part = GetVtk();
602  vtkSmartPointer<vtkPoints> p_points = vtkSmartPointer<vtkPoints>::New();
603 
604  for(unsigned idx=0; idx<rLocations.size(); idx++)
605  {
606  if(DIM==3)
607  {
608  p_points->InsertNextPoint(rLocations[idx].GetLocation(mReferenceLength)[0],
609  rLocations[idx].GetLocation(mReferenceLength)[1], rLocations[idx].GetLocation(mReferenceLength)[2]);
610  }
611  else
612  {
613  p_points->InsertNextPoint(rLocations[idx].GetLocation(mReferenceLength)[0],
614  rLocations[idx].GetLocation(mReferenceLength)[1], 0.0);
615  }
616  }
617 
618  vtkSmartPointer<vtkPolyData> p_point_data = vtkSmartPointer<vtkPolyData>::New();
619  p_point_data->SetPoints(p_points);
620 
621  //Points inside test
622  vtkSmartPointer<vtkSelectEnclosedPoints> selectEnclosedPoints = vtkSmartPointer<vtkSelectEnclosedPoints>::New();
623  #if VTK_MAJOR_VERSION <= 5
624  selectEnclosedPoints->SetInput(p_point_data);
625  #else
626  selectEnclosedPoints->SetInputData(p_point_data);
627  #endif
628  #if VTK_MAJOR_VERSION <= 5
629  selectEnclosedPoints->SetSurface(p_part);
630  #else
631  selectEnclosedPoints->SetSurfaceData(p_part);
632  #endif
633  selectEnclosedPoints->Update();
634 
635  std::vector<bool> is_inside(rLocations.size());
636  for(unsigned idx=0; idx<rLocations.size(); idx++)
637  {
638  is_inside[idx] = selectEnclosedPoints->IsInside(idx);
639  }
640  return is_inside;
641 }
642 
643 template<unsigned DIM>
645 {
646  // Loop through the nodes of each polygon. If it is in another polygon, replace it.
647  std::vector<boost::shared_ptr<Polygon<DIM> > > polygons = GetPolygons();
648  for(unsigned idx=0; idx<polygons.size(); idx++)
649  {
650  for(unsigned jdx=0; jdx<polygons.size(); jdx++)
651  {
652  if(idx != jdx)
653  {
654  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > p1_verts = polygons[idx]->GetVertices();
655  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > p2_verts = polygons[jdx]->GetVertices();
656  for(unsigned mdx=0; mdx<p1_verts.size(); mdx++)
657  {
658  for(unsigned ndx=0; ndx<p2_verts.size(); ndx++)
659  {
660  if(GetDistance(*(p2_verts[ndx]), *(p1_verts[mdx]))< 1.e-6*p2_verts[ndx]->GetReferenceLengthScale())
661  {
662  polygons[jdx]->ReplaceVertex(ndx, polygons[idx]->GetVertex(mdx));
663  }
664  }
665  }
666  }
667  }
668  }
669 
670  for(unsigned idx=0; idx<mFacets.size(); idx++)
671  {
672  mFacets[idx]->UpdateVertices();
673  }
674  mVtkIsUpToDate = false;
675 }
676 
677 template<unsigned DIM>
678 void Part<DIM>::SetReferenceLengthScale(units::quantity<unit::length> referenceLength)
679 {
680  mReferenceLength = referenceLength;
681  mVtkIsUpToDate = false;
682 }
683 
684 template<unsigned DIM>
686 {
687  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > vertices = GetVertices();
688  {
689  for (unsigned idx = 0; idx < vertices.size(); idx++)
690  {
691  vertices[idx]->Translate(vector);
692  }
693  }
694  mVtkIsUpToDate = false;
695 }
696 
697 template<unsigned DIM>
698 void Part<DIM>::Write(const std::string& fileName, GeometryFormat::Value format)
699 {
700  GeometryWriter writer;
701  writer.SetFileName(fileName);
702  writer.SetInput(GetVtk());
703  writer.SetOutputFormat(format);
704  writer.Write();
705 }
706 
707 // Explicit instantiation
708 template class Part<2> ;
709 template class Part<3> ;
void SetReferenceLengthScale(units::quantity< unit::length > referenceLength)
Definition: Part.cpp:678
std::vector< unsigned > GetContainingGridIndices(unsigned num_x, unsigned num_y, unsigned num_z, units::quantity< unit::length > spacing)
Definition: Part.cpp:386
void MergeCoincidentVertices()
Definition: Part.cpp:644
vtkSmartPointer< vtkPolyData > GetVtk()
Definition: Part.cpp:511
void BooleanWithNetwork(boost::shared_ptr< VesselNetwork< DIM > > pVesselNetwork)
Definition: Part.cpp:323
units::quantity< unit::length > GetReferenceLengthScale() const
boost::shared_ptr< Facet< DIM > > GetFacet(const DimensionalChastePoint< DIM > &rLocation)
Definition: Part.cpp:351
bool IsPointInPart(DimensionalChastePoint< DIM > location)
Definition: Part.cpp:564
void AddCylinder(units::quantity< unit::length > radius, units::quantity< unit::length > depth, DimensionalChastePoint< DIM > centre, unsigned numSegments=24)
Definition: Part.cpp:100
boost::shared_ptr< Polygon< DIM > > AddRectangle(units::quantity< unit::length > sizeX, units::quantity< unit::length > sizeY, DimensionalChastePoint< DIM > origin)
Definition: Part.cpp:161
std::vector< DimensionalChastePoint< DIM > > GetVertexLocations()
Definition: Part.cpp:409
void SetOutputFormat(GeometryFormat::Value format)
Part()
Definition: Part.cpp:52
std::vector< std::pair< unsigned, unsigned > > GetSegmentIndices()
Definition: Part.cpp:478
std::vector< DimensionalChastePoint< DIM > > GetHoleMarkers()
Definition: Part.cpp:339
std::vector< boost::shared_ptr< DimensionalChastePoint< DIM > > > GetVertices()
Definition: Part.cpp:364
std::vector< boost::shared_ptr< Facet< DIM > > > GetFacets()
Definition: Part.cpp:472
void Write(const std::string &rFilename, GeometryFormat::Value format=GeometryFormat::VTP)
Definition: Part.cpp:698
units::quantity< unit::length > mReferenceLength
Definition: Part.hpp:85
void Extrude(boost::shared_ptr< Polygon< DIM > > pPolygon, units::quantity< unit::length > distance)
Definition: Part.cpp:274
void SetInput(vtkSmartPointer< vtkPolyData > pSurface)
void AddCuboid(units::quantity< unit::length > sizeX, units::quantity< unit::length > sizeY, units::quantity< unit::length > sizeZ, DimensionalChastePoint< DIM > origin)
Definition: Part.cpp:111
std::vector< boost::shared_ptr< Facet< DIM > > > mFacets
Definition: Part.hpp:65
std::vector< DimensionalChastePoint< DIM > > mHoleMarkers
Definition: Part.hpp:75
void Translate(DimensionalChastePoint< DIM > vector)
Definition: Part.cpp:685
static boost::shared_ptr< Polygon > Create(std::vector< boost::shared_ptr< DimensionalChastePoint< DIM > > > vertices)
Definition: Polygon.cpp:60
std::vector< boost::shared_ptr< Polygon< DIM > > > GetSurfacePolygons()
boost::shared_ptr< Polygon< DIM > > AddPolygon(std::vector< boost::shared_ptr< DimensionalChastePoint< DIM > > > vertices, bool newFacet=false, boost::shared_ptr< Facet< DIM > > pFacet=boost::shared_ptr< Facet< DIM > >())
Definition: Part.cpp:128
units::quantity< unit::length > GetReferenceLengthScale()
Definition: Part.cpp:345
bool mVtkIsUpToDate
Definition: Part.hpp:90
c_vector< double, DIM > GetLocation(units::quantity< unit::length > scale)
std::vector< DimensionalChastePoint< DIM > > GetHoles()
boost::shared_ptr< Polygon< DIM > > AddCircle(units::quantity< unit::length > radius, DimensionalChastePoint< DIM > centre, unsigned numSegments=24)
Definition: Part.cpp:76
Definition: Facet.hpp:54
vtkSmartPointer< vtkPolyData > mVtkPart
Definition: Part.hpp:70
std::vector< boost::shared_ptr< Polygon< DIM > > > GetPolygons()
Definition: Part.cpp:421
~Part()
Definition: Part.cpp:70
std::vector< units::quantity< unit::length > > GetBoundingBox()
Definition: Part.cpp:433
static boost::shared_ptr< Part< DIM > > Create()
Definition: Part.cpp:63
void AddVesselNetwork(boost::shared_ptr< VesselNetwork< DIM > > pVesselNetwork, bool surface=false)
Definition: Part.cpp:189
void SetFileName(const std::string &rFileName)
void AddHoleMarker(DimensionalChastePoint< DIM > location)
Definition: Part.cpp:122
Definition: Part.hpp:60