intersect_polys_r3d.h
Go to the documentation of this file.
1 /*
2 This file is part of the Ristra portage project.
3 Please see the license file at the root of this repository, or at:
4  https://github.com/laristra/portage/blob/master/LICENSE
5 */
6 
7 #ifndef INTERSECT_POLYS_R3D_H
8 #define INTERSECT_POLYS_R3D_H
9 
10 #include <array>
11 #include <stdexcept>
12 #include <vector>
13 #include <algorithm>
14 
15 extern "C" {
16 #include "wonton/intersect/r3d/r3d.h"
17 }
18 
19 #ifdef HAVE_TANGRAM
20 #include "tangram/support/MatPoly.h"
21 #endif
22 
23 
25 #include "wonton/support/Point.h"
26 #include "wonton/support/Vector.h"
27 
28 
29 namespace Portage {
30 
31 using Wonton::Point;
32 using Wonton::Vector;
33 
34 typedef
35 struct facetedpoly {
36  std::vector<std::vector<int>> facetpoints; // we can flatten this to a
37  // // 1D array if we assume only
38  // // triangular facets or we
39  // // include the number of points
40  // // in each facet
41  std::vector<Point<3>> points;
43 
44 #ifdef HAVE_TANGRAM
45 
50 inline
51 facetedpoly_t get_faceted_matpoly(const Tangram::MatPoly<3>& matpoly) {
52  // facet the matpoly
53  Tangram::MatPoly<3> faceted_matpoly;
54  matpoly.faceted_matpoly(&faceted_matpoly);
55 
56  // initialize the result with the face vertices
57  facetedpoly_t result{faceted_matpoly.face_vertices()};
58 
59  // convert tangram points to portage points for facetedpoly_t.points
60  for (auto p : faceted_matpoly.points()) result.points.push_back(Point<3>(p));
61  return result;
62 }
63 #endif
64 
65 
66 // Intersect one source polyhedron (possibly non-convex but with
67 // triangular facets only) with a bunch of tets forming a target
68 // polyhedron
69 
70 std::vector<double>
71 inline
73  const std::vector<std::array<Point<3>, 4>> &target_tet_coords,
74  NumericTolerances_t num_tols) {
75 
76  // Bounding box of the target cell - will be used to compute
77  // epsilon for bounding box check. We could use the source cell
78  // coordinates to find the bounds as well but its probably an
79  // unnecessary computation (EVEN if the target cell is much
80  // smaller than the source cell)
81 
82  std::vector<double> source_cell_bounds = {1e99, -1e99, 1e99, -1e99,
83  1e99, -1e99};
84 
85  // Initialize the source polyhedron description in a form R3D wants
86  // Simultaneously compute the bounding box
87  r3d_poly src_r3dpoly;
88  int num_verts = srcpoly.points.size();
89  r3d_rvec3 *verts = new r3d_rvec3[num_verts];
90  for (int i = 0; i < num_verts; i++) {
91  for (int j = 0; j < 3; j++) {
92  verts[i].xyz[j] = srcpoly.points[i][j];
93 
94  if (source_cell_bounds[2*j] > verts[i].xyz[j])
95  source_cell_bounds[2*j] = verts[i].xyz[j];
96  if (source_cell_bounds[2*j+1] < verts[i].xyz[j])
97  source_cell_bounds[2*j+1] = verts[i].xyz[j];
98  }
99  }
100 
101  // used only for bounding box check not for intersections
102  double bbeps = num_tols.min_absolute_distance;
103 
104  int num_faces = srcpoly.facetpoints.size();
105  r3d_int *face_num_verts = new r3d_int[num_faces];
106  for (int i = 0; i < num_faces; i++)
107  face_num_verts[i] = srcpoly.facetpoints[i].size();
108 
109  r3d_int **face_vert_ids = new r3d_int *[num_faces];
110  for (int i = 0; i < num_faces; i++)
111  face_vert_ids[i] = new r3d_int[face_num_verts[i]];
112 
113  for (int i = 0; i < num_faces; i++)
114  std::copy(srcpoly.facetpoints[i].begin(), srcpoly.facetpoints[i].end(),
115  face_vert_ids[i]);
116 
117 #ifdef DEBUG
118  // Lets check the volume of the source polygon - If its convex or
119  // mildly non-convex, we should get a positive volume
120 
121  // First calculate a center point
122  Point<3> cen(0.0, 0.0, 0.0);
123  for (int i = 0; i < num_verts; i++)
124  cen += srcpoly.points[i];
125  cen /= num_verts;
126 
127  // Assume triangular facets only
128  double polyvol = 0.0;
129  for (int i = 0; i < num_faces; i++) {
130  // p0, p1, p2 traversed in order form a triangle whose normal
131  // points out of the source polyhedron
132  const Point<3> &p0 = srcpoly.points[srcpoly.facetpoints[i][0]];
133  const Point<3> &p1 = srcpoly.points[srcpoly.facetpoints[i][1]];
134  const Point<3> &p2 = srcpoly.points[srcpoly.facetpoints[i][2]];
135 
136  Vector<3> v0 = p1-p0;
137  Vector<3> v1 = p2-p0;
138  Vector<3> outnormal = cross(v0, v1);
139  Vector<3> v2 = cen-p0;
140  double tetvol = -dot(v2, outnormal)/6.0;
141  if (tetvol < 0.0)
142  std::cerr << "Wrong orientation of facets or non-convex polyhedron" <<
143  std::endl;
144  polyvol += tetvol;
145  }
146  if (polyvol < 0.0)
147  throw std::runtime_error("Source polyhedron has negative volume");
148 #endif
149 
150  r3d_init_poly(&src_r3dpoly, verts, num_verts, face_vert_ids, face_num_verts,
151  num_faces);
152 
153  // Finished building source poly; now intersect with tets of target cell
154 
155  std::vector<double> moments(4, 0);
156  for (auto const & target_cell_tet : target_tet_coords) {
157  std::vector<r3d_plane> faces(4);
158 
159  std::vector<double> target_tet_bounds = {1e99, -1e99, 1e99,
160  -1e99, 1e99, -1e99};
161  r3d_rvec3 verts2[4];
162  for (int i = 0; i < 4; i++)
163  for (int j = 0; j < 3; j++) {
164  verts2[i].xyz[j] = target_cell_tet[i][j];
165 
166  if (target_tet_bounds[2*j] > verts2[i].xyz[j])
167  target_tet_bounds[2*j] = verts2[i].xyz[j];
168  if (target_tet_bounds[2*j+1] < verts2[i].xyz[j])
169  target_tet_bounds[2*j+1] = verts2[i].xyz[j];
170  }
171 
172  // Check if the target and source bounding boxes overlap - bbeps
173  // is used to subject touching tets to the full intersection
174  // (just in case)
175 
176  bool check_bb = true;
177  if (check_bb) {
178  bool disjoint = false;
179  for (int j = 0; j < 3 && !disjoint; ++j)
180  disjoint = (target_tet_bounds[2*j] > source_cell_bounds[2*j+1]+bbeps ||
181  target_tet_bounds[2*j+1] < source_cell_bounds[2*j]-bbeps);
182  if (disjoint) continue;
183  }
184 
185 #ifdef DEBUG
186  if (r3d_orient(verts2) < 0)
187  throw std::runtime_error("target_wedge has negative volume");
188 #endif
189 
190  r3d_tet_faces_from_verts(&faces[0], verts2);
191 
192  // clip the source poly against the faces of the target tet - but
193  // make a copy of src_r3dpoly first because it will get modified
194  // in the process of clipping
195  r3d_poly src_r3dpoly_copy = src_r3dpoly;
196  r3d_clip(&src_r3dpoly_copy, &faces[0], 4);
197 
198  // find the moments (up to quadratic order) of the clipped poly
199  const int POLY_ORDER = 1;
200  r3d_real om[R3D_NUM_MOMENTS(POLY_ORDER)];
201  r3d_reduce(&src_r3dpoly_copy, om, POLY_ORDER);
202 
203  // Check that the returned volume is positive (if the volume is
204  // zero, i.e. abs(om[0]) < eps, then it can sometimes be
205  // slightly negative, like om[0] == -1.24811e-16.
206  // @todo multiply by domain or element size
207  if (om[0] < num_tols.minimal_intersection_volume)
208  throw std::runtime_error("Negative volume");
209 
210  // Accumulate moments:
211  for (int i = 0; i < 4; i++)
212  moments[i] += om[i];
213  }
214 
215  delete [] verts;
216  delete [] face_num_verts;
217  for (int i = 0; i < num_faces; i++)
218  delete [] face_vert_ids[i];
219  delete [] face_vert_ids;
220 
221  return moments;
222 } // intersect_polys_3D
223 
224 } // namespace Portage
225 
226 #endif // INTERSECT_POLYS_R3D_H
double min_absolute_distance
Definition: portage.h:165
std::vector< T > vector
Definition: portage.h:238
Intersection and other tolerances to handle tiny values.
Definition: portage.h:152
double minimal_intersection_volume
Definition: portage.h:160
Definition: intersect_polys_r3d.h:34
std::vector< double > intersect_polys_r3d(const facetedpoly_t &srcpoly, const std::vector< std::array< Point< 3 >, 4 >> &target_tet_coords, NumericTolerances_t num_tols)
Definition: intersect_polys_r3d.h:72
Definition: coredriver.h:42
std::vector< Point< 3 > > points
Definition: intersect_polys_r3d.h:41
std::vector< std::vector< int > > facetpoints
Definition: intersect_polys_r3d.h:36
struct Portage::facetedpoly facetedpoly_t