b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2overlaying_surface_meshes.H
1//------------------------------------------------------------------------
2// b2overlaying_surface_meshes.H --
3//
4// written by doreille
5//
6// Copyright (c) 2010-2012 SMR Engineering & Development SA
7// 2502 Bienne, Switzerland
8//
9// All Rights Reserved. Proprietary source code. The contents of
10// this file may not be disclosed to third parties, copied or
11// duplicated in any form, in whole or in part, without the prior
12// written permission of SMR.
13//------------------------------------------------------------------------
14
15#ifndef B2OVERLAYING_SURFACE_MESHES_H_
16#define B2OVERLAYING_SURFACE_MESHES_H_
17
18#include <algorithm>
19#include <cassert>
20#include <deque>
21#include <limits>
22#include <list>
23#include <stack>
24#include <vector>
25
26#include "b2exception.H"
27#include "b2nonlinear_system_solver.H"
28#include "b2tensor_calculus.H"
29
30// #define B2OSM_DEBUG 1
31// #define AVOID_PROBLEM_IN_TRIANGLE_ITERATOR 1
32
33namespace b2000 {
34
35/* prototype class for the template compute_overlaying_surface_meshes function
36
37 struct B_Vertex {
38 DCELHalfEdge* edge;
39 double coor[3];
40 double internal_coor[2];
41 void* subvertex;
42 };
43
44 struct GVertex {
45 DCELHalfEdge* edge;
46 double coor[3];
47 double internal_coor[2];
48 double normal[3];
49 void* subvertex;
50 };
51
52 struct Face {
53 DCELHalfEdge* edge;
54 };
55
56 struct B_DCELHalfEdge {
57 B_Vertex* vertex; // vertex at the end of the half-edge
58 DCELHalfEdge* pair; // oppositely oriented adjacent half-edge
59 Face* face; // face the half-edge borders
60 DCELHalfEdge* next; // next half-edge around the face
61 void* subvertex;
62 };
63
64 struct G_DCELHalfEdge {
65 G_Vertex* vertex; // vertex at the end of the half-edge
66 G_DCELHalfEdge* pair; // oppositely oriented adjacent half-edge
67 Face* face; // face the half-edge borders
68 G_DCELHalfEdge* next; // next half-edge around the face
69 void* subvertex;
70 };
71
72 // the border edges must also have a pair edge.
73 struct DCELMesh {
74 typedef vertex_type;
75 typedef vertex_iterator;
76
77 typedef face_type;
78 typedef face_iterator;
79
80 typedef edge_type;
81 typedef edge_iterator;
82
83 vertex_const_iterator vertex_begin();
84 vertex_const_iterator vertex_end();
85 size_t vertex_size() const;
86
87 face_const_iterator face_begin();
88 face_const_iterator face_end();
89 size_t face_size() const;
90
91 edge_const_iterator edge_begin();
92 edge_const_iterator edge_end();
93 size_t edge_size() const;
94 };
95
96 struct Triangle {
97 Face* face_b;
98 double internal_coor_b[3][2];
99 Face* face_g;
100 double internal_coor_g[3][2];
101 };
102
103 */
104
112//
113// The algorithm is described in the paper:
114// Overlaying Surface Meshes, Part I: Algorithms,
115// Xiangmin Jiao and Michael T. Heath,
116// International Journal of Computational Geometry & Applications,
117// Vol. 14, No 6 (2004) 379-402
119// Internal class that compute the overlaying surface meshes
120template <typename G_DCEL_MESH, typename B_DCEL_MESH>
121class ComputeOverlayingSurfaceMeshes {
122private:
123 struct SubVertex;
124
125public:
126 typedef typename G_DCEL_MESH::edge_type g_edge_type;
127 typedef typename G_DCEL_MESH::vertex_type g_vertex_type;
128 typedef typename G_DCEL_MESH::face_type g_face_type;
129 typedef typename B_DCEL_MESH::edge_type b_edge_type;
130 typedef typename B_DCEL_MESH::vertex_type b_vertex_type;
131 typedef typename B_DCEL_MESH::face_type b_face_type;
132
133 void init(G_DCEL_MESH& g_mesh_, B_DCEL_MESH& b_mesh_, double tolerancing_ = 1e-5) {
134 g_mesh = &g_mesh_;
135 b_mesh = &b_mesh_;
136 tolerancing = tolerancing_;
137 epsilon = 1e-12;
138 g_edge_type* g_border_edge;
139
140 // initialisation of subvertex of green and blue edge to NULL
141 g_border_edge = 0;
142 for (typename G_DCEL_MESH::edge_iterator i = g_mesh->edge_begin(); i != g_mesh->edge_end();
143 ++i) {
144 i->subvertex = 0;
145 if (i->face == 0) { g_border_edge = &*i; }
146 }
147 b_border_edge = 0;
148 for (typename B_DCEL_MESH::edge_iterator i = b_mesh->edge_begin(); i != b_mesh->edge_end();
149 ++i) {
150 i->subvertex = 0;
151 if (i->face == 0) { b_border_edge = &*i; }
152 }
153
154 // step one: locating subvertices along blue edges
155 {
156 std::stack<b_edge_type*> edge_to_process;
157
158 // brute force to find a first green parent of a blue vertex
159 {
160 std::pair<b_vertex_type*, GSubVertex> fsv =
161 green_parent_of_a_blue_vertex_by_brute_force();
162 SubVertex& sv = new_sub_vertex();
163 sv.g_subvertex = fsv.second;
164 sv.b_subvertex.dim = 0;
165 sv.b_subvertex.v = fsv.first;
166 fsv.first->subvertex = &sv;
167
168 b_vertex_iterator_on_edge iter(fsv.first);
169 for (b_edge_type* e = iter.next(); e != 0; e = iter.next()) {
170 edge_to_process.push(e);
171 }
172 }
173
174 // iteration on the blue edge in a breadth-first order
175 std::vector<g_edge_type*> edge_candidate;
176 edge_candidate.reserve(20);
177 std::vector<g_edge_type*> edge_candidate_tmp;
178 edge_candidate_tmp.reserve(20);
179 while (!edge_to_process.empty()) {
180 b_edge_type* b_edge = edge_to_process.top();
181 // std::cout << "b edge to process = " << b_edge->pair->vertex->id + 1 << "->" <<
182 // b_edge->vertex->id + 1 << std::endl;
183 SubVertex* g_i1 = 0;
184 SubVertex* g_i2 = static_cast<SubVertex*>(b_edge->pair->vertex->subvertex);
185 SubVertex* g_i3;
186 edge_to_process.pop();
187 if (b_edge->subvertex) { continue; }
188 for (;;) {
189 edge_candidate.clear();
190 switch (g_i2->g_subvertex.dim) {
191 case 0: {
192 g_edge_type* e_start = g_i2->g_subvertex.v->edge->next;
193 g_edge_type* e = e_start;
194 do {
195 bool in_g_i1 = false;
196 if (g_i1 != 0) {
197 switch (g_i1->g_subvertex.dim) {
198 case 0: {
199 g_vertex_type* v = g_i1->g_subvertex.v;
200 in_g_i1 = v == e->vertex || v == e->pair->vertex;
201 break;
202 }
203 case 1:
204 in_g_i1 = g_i1->g_subvertex.e == e
205 || g_i1->g_subvertex.e == e->pair;
206 break;
207 case 2:
208 in_g_i1 = e->face == g_i1->g_subvertex.f;
209 }
210 }
211 bool ext_but_on_g_i2_face = false;
212 if (e->face == 0) {
213 g_vertex_type* v = g_i2->g_subvertex.v;
214 g_edge_type* e_tmp = e->pair;
215 do {
216 if (e_tmp->vertex == v) {
217 ext_but_on_g_i2_face = true;
218 break;
219 }
220 e_tmp = e_tmp->next;
221 } while (e_tmp != e->pair);
222 }
223 if (!in_g_i1 && !ext_but_on_g_i2_face) {
224 edge_candidate.push_back(e);
225 }
226 e = e->next;
227 if (e->vertex == g_i2->g_subvertex.v) { e = e->pair->next; }
228 } while (e != e_start);
229 break;
230 }
231 case 1: {
232 g_edge_type* e_start = g_i2->g_subvertex.e->pair;
233 if (g_i1 == 0 && e_start->face == 0) { e_start = e_start->pair; }
234 g_edge_type* e = e_start->next;
235 do {
236 edge_candidate.push_back(e);
237 e = e->next;
238 } while (e != e_start);
239 if (g_i1 == 0) {
240 e_start = e_start->pair;
241 e = e_start->next;
242 do {
243 if (e->pair->face != e_start->pair->face) {
244 edge_candidate.push_back(e);
245 }
246 e = e->next;
247 } while (e != e_start);
248 }
249 break;
250 }
251 case 2: {
252 g_edge_type* e_start =
253 g_i2->g_subvertex.f ? g_i2->g_subvertex.f->edge : g_border_edge;
254 g_edge_type* e = e_start;
255 do {
256 edge_candidate.push_back(e);
257 e = e->next;
258 } while (e != e_start);
259 break;
260 }
261 }
262
263 // reorder the edge candidates list to put the external edge at the end of the
264 // list
265 {
266 size_t i = 0;
267 while (i != edge_candidate.size() && edge_candidate[i]->face != 0) { ++i; }
268 while (i != edge_candidate.size() && edge_candidate[i]->face == 0) { ++i; }
269 if (i != edge_candidate.size()) {
270 edge_candidate_tmp.assign(
271 edge_candidate.begin() + i, edge_candidate.end());
272 edge_candidate_tmp.insert(
273 edge_candidate_tmp.end(), edge_candidate.begin(),
274 edge_candidate.begin() + i);
275 edge_candidate.swap(edge_candidate_tmp);
276 }
277 }
278
279 g_i3 = 0;
280 bool first_res = false;
281 double first_alpha = 0;
282 double first_beta = 0;
283 SubVertex first_sv;
284 double last_alpha = -1;
285 double last_beta = -1;
286 bool last_res = false;
287 SubVertex sv;
288 SubVertex last_sv;
289 SubVertex best_sv;
290 SubVertex* best_g_i3 = 0;
291 double best_alpha = 2;
292 size_t i_end = edge_candidate.size();
293 if (edge_candidate.size() > 1
294 && edge_candidate.front()->pair->vertex == edge_candidate.back()->vertex) {
295 ++i_end;
296 }
297 for (size_t i = 0; i != i_end; ++i) {
298 double alpha = 0;
299 double beta = 0;
300 bool res;
301 if (i == edge_candidate.size()) {
302 res = first_res;
303 sv = first_sv;
304 alpha = first_alpha;
305 beta = first_beta;
306 } else {
307 res = edge_intersection(*b_edge, *edge_candidate[i], sv, alpha, beta);
308 if (res && g_i2 && g_i2->b_subvertex.dim == 1
309 && alpha < g_i2->b_subvertex.icoor[0]) {
310 res = false;
311 }
312 }
313 if (i == 0) {
314 first_res = res;
315 first_sv = sv;
316 first_alpha = alpha;
317 first_beta = beta;
318 }
319
320 if (last_res) {
321 if (res) {
322 if (std::abs(last_alpha - alpha) < tolerancing) {
323 // case (b) of section 6.1 of the reference article
324 if (sv.b_subvertex.dim == 0) {
325 if (sv.b_subvertex.v->subvertex == 0) {
326 g_i3 = &sv;
327 sv.b_subvertex.v = b_edge->vertex;
328 } else {
329 g_i3 = static_cast<SubVertex*>(
330 sv.b_subvertex.v->subvertex);
331 }
332 } else if (last_sv.b_subvertex.dim == 0) {
333 if (last_sv.b_subvertex.v->subvertex == 0) {
334 g_i3 = &sv;
335 sv.b_subvertex.dim = 0;
336 sv.b_subvertex.v = b_edge->vertex;
337 } else {
338 g_i3 = static_cast<SubVertex*>(
339 last_sv.b_subvertex.v->subvertex);
340 }
341 } else {
342 g_i3 = &sv;
343 sv.b_subvertex.dim = 1;
344 sv.b_subvertex.e = b_edge;
345 alpha = (last_alpha + alpha) / 2;
346 sv.b_subvertex.icoor[0] = alpha;
347 }
348 g_i3->g_subvertex.dim = 0;
349 g_i3->g_subvertex.v =
350 edge_candidate[i - (last_beta > beta ? 1 : 0)]->vertex;
351 } else {
352 if (alpha < last_alpha) {
353 if (sv.b_subvertex.dim == 0
354 && sv.b_subvertex.v->subvertex) {
355 // case (a) of section 6.1 of the reference article
356 g_i3 = static_cast<SubVertex*>(
357 sv.b_subvertex.v->subvertex);
358 if (g_i3->g_subvertex.dim >= last_sv.g_subvertex.dim) {
359 if (g_i3->g_subvertex.dim == 1
360 && last_sv.g_subvertex.dim == 1) {
361 g_i3->g_subvertex.dim = 0;
362 if (g_i3->g_subvertex.icoor[0] < 0.5) {
363 g_i3->g_subvertex.v =
364 g_i3->g_subvertex.e->pair->vertex;
365 } else {
366 g_i3->g_subvertex.v =
367 g_i3->g_subvertex.e->vertex;
368 }
369 } else {
370 g_i3->g_subvertex = sv.g_subvertex;
371 }
372 }
373 } else {
374 g_i3 = &sv;
375 }
376 } else {
377 if (last_sv.b_subvertex.dim == 0
378 && last_sv.b_subvertex.v->subvertex) {
379 // case (a) of section 6.1 of the reference article
380 g_i3 = static_cast<SubVertex*>(
381 last_sv.b_subvertex.v->subvertex);
382 if (g_i3->g_subvertex.dim >= last_sv.g_subvertex.dim) {
383 if (g_i3->g_subvertex.dim == 1
384 && last_sv.g_subvertex.dim == 1) {
385 g_i3->g_subvertex.dim = 0;
386 if (g_i3->g_subvertex.icoor[0] < 0.5) {
387 g_i3->g_subvertex.v =
388 g_i3->g_subvertex.e->pair->vertex;
389 } else {
390 g_i3->g_subvertex.v =
391 g_i3->g_subvertex.e->vertex;
392 }
393 } else {
394 g_i3->g_subvertex = last_sv.g_subvertex;
395 }
396 }
397 } else {
398 g_i3 = &sv;
399 sv.g_subvertex = last_sv.g_subvertex;
400 sv.b_subvertex = last_sv.b_subvertex;
401 }
402 }
403 }
404 } else {
405 if (last_sv.b_subvertex.dim == 0
406 && last_sv.b_subvertex.v->subvertex) {
407 // case (a) of section 6.1 of the reference article
408 g_i3 =
409 static_cast<SubVertex*>(last_sv.b_subvertex.v->subvertex);
410 if (g_i3->g_subvertex.dim >= last_sv.g_subvertex.dim) {
411 g_i3->g_subvertex = last_sv.g_subvertex;
412 }
413 } else {
414 g_i3 = &sv;
415 sv.g_subvertex = last_sv.g_subvertex;
416 sv.b_subvertex = last_sv.b_subvertex;
417 }
418 }
419 } else if (!res && last_beta > 1 && beta < 1) {
420 double max_alpha = std::max(last_alpha, alpha);
421 double min_alpha = std::min(last_alpha, alpha);
422 alpha = (last_alpha + alpha) / 2;
423 if (min_alpha > 1 - 2 * tolerancing && max_alpha < 1 + tolerancing) {
424 // case (c) of section 6.1 of the reference article
425 if (max_alpha < tolerancing) {
426 g_i3 = static_cast<SubVertex*>(b_edge->pair->vertex->subvertex);
427 } else if (min_alpha < 1 - tolerancing) {
428 g_i3 = &sv;
429 sv.b_subvertex.dim = 0;
430 sv.b_subvertex.v = b_edge->vertex;
431 } else {
432 if (b_edge->vertex->subvertex) {
433 g_i3 = static_cast<SubVertex*>(b_edge->vertex->subvertex);
434 } else {
435 g_i3 = &sv;
436 g_i3->b_subvertex.dim = 1;
437 g_i3->b_subvertex.e = b_edge;
438 g_i3->b_subvertex.icoor[0] = alpha;
439 }
440 }
441 g_i3->g_subvertex.dim = 0;
442 g_i3->g_subvertex.v = edge_candidate[i - 1]->vertex;
443 }
444 }
445 if (g_i2->g_subvertex.f == 0) {
446 if (g_i3 && g_i3->b_subvertex.dim == 1
447 && g_i3->b_subvertex.icoor[0] < best_alpha) {
448 if (g_i3 == &sv) {
449 best_sv = sv;
450 best_g_i3 = &best_sv;
451 } else {
452 best_g_i3 = g_i3;
453 }
454 best_alpha = g_i3->b_subvertex.icoor[0];
455 }
456 g_i3 = 0;
457 } else if (g_i3) {
458 break;
459 }
460 last_alpha = alpha;
461 last_beta = beta;
462 last_res = res;
463 last_sv = sv;
464 }
465
466 if (g_i3 == 0 && last_res && edge_candidate.size() == i_end) {
467 if (last_sv.b_subvertex.dim == 0 && last_sv.b_subvertex.v->subvertex) {
468 // case (a) of section 6.1 of the reference article
469 g_i3 = static_cast<SubVertex*>(last_sv.b_subvertex.v->subvertex);
470 if (g_i3->g_subvertex.dim >= last_sv.g_subvertex.dim) {
471 g_i3->g_subvertex = last_sv.g_subvertex;
472 }
473 } else {
474 g_i3 = &sv;
475 sv.g_subvertex = last_sv.g_subvertex;
476 sv.b_subvertex = last_sv.b_subvertex;
477 }
478 if (g_i2->g_subvertex.f == 0 && g_i3->b_subvertex.icoor[0] < best_alpha) {
479 if (g_i3 == &sv) {
480 best_sv = sv;
481 best_g_i3 = &best_sv;
482 } else {
483 best_g_i3 = g_i3;
484 }
485 best_alpha = g_i3->b_subvertex.icoor[0];
486 }
487 }
488
489 if (g_i2->g_subvertex.f == 0) {
490 if (best_g_i3 == 0) {
491 if (b_edge->vertex->subvertex == 0) {
492 SubVertex& rsv = new_sub_vertex();
493 rsv.b_subvertex.dim = 0;
494 rsv.b_subvertex.v = b_edge->vertex;
495 rsv.g_subvertex.dim = 2;
496 rsv.g_subvertex.f = 0;
497 std::fill_n(rsv.g_subvertex.icoor, 2, 0.0);
498 g_i3 = &rsv;
499 } else {
500 g_i3 = static_cast<SubVertex*>(b_edge->vertex->subvertex);
501 }
502 } else {
503 if (best_g_i3 == &best_sv) {
504 SubVertex& sv = new_sub_vertex();
505 sv.g_subvertex = best_g_i3->g_subvertex;
506 sv.b_subvertex = best_g_i3->b_subvertex;
507 g_i3 = &sv;
508 } else {
509 g_i3 = best_g_i3;
510 }
511 }
512 } else {
513 if (g_i3 == &sv) {
514 SubVertex& sv = new_sub_vertex();
515 sv.g_subvertex = g_i3->g_subvertex;
516 sv.b_subvertex = g_i3->b_subvertex;
517 g_i3 = &sv;
518 }
519 if (g_i3 == 0) {
520 // the blue edge do not intercept a green edge
521 if (b_edge->vertex->subvertex == 0) {
522 SubVertex& rsv = new_sub_vertex();
523 rsv.b_subvertex.dim = 0;
524 rsv.b_subvertex.v = b_edge->vertex;
525 rsv.b_subvertex.icoor[0] = rsv.b_subvertex.icoor[1] = 0;
526 std::set<g_face_type*> processed_face;
527 for (size_t i = 0; i != edge_candidate.size(); ++i) {
528 g_face_type* f = edge_candidate[i]->face;
529 if (f == 0) { continue; }
530 if (processed_face.insert(f).second) {
531 double dist;
532 if (b_point_to_g_face(
533 *rsv.b_subvertex.v, *f, rsv.g_subvertex, dist)) {
534 g_i3 = &rsv;
535 break;
536 }
537 }
538 }
539 if (g_i3 == 0) {
540 // the blue vertex do not project on the green mesh
541 rsv.g_subvertex.dim = 2;
542 rsv.g_subvertex.f = 0;
543 std::fill_n(rsv.g_subvertex.icoor, 2, 0.0);
544 g_i3 = &rsv;
545 }
546 } else {
547 g_i3 = static_cast<SubVertex*>(b_edge->vertex->subvertex);
548 }
549 }
550 }
551
552 const int g_i2_b_subvertex_dim = g_i2->b_subvertex.dim;
553 if (g_i3->g_subvertex.dim == 0) {
554 if (g_i3->g_subvertex.v->subvertex == 0) {
555 g_i3->g_subvertex.v->subvertex = g_i3;
556 } else if (g_i3 == &list_subvertex.back()) {
557 SubVertex* tmp =
558 static_cast<SubVertex*>(g_i3->g_subvertex.v->subvertex);
559 if (tmp->b_subvertex.dim != 0 || g_i3->b_subvertex.dim != 0) {
560 cluster_subvertex_on_dcel(tmp, g_i3);
561 list_subvertex.pop_back();
562 g_i3 = tmp;
563 }
564 }
565 }
566 if (g_i3 != g_i2) {
567 if (b_edge < b_edge->pair) {
568 g_i3->b_previous_subvertex = g_i2;
569 g_i2->b_next_subvertex = g_i3;
570 } else {
571 g_i3->b_next_subvertex = g_i2;
572 g_i2->b_previous_subvertex = g_i3;
573 }
574 }
575 if (g_i3->b_subvertex.dim == 0) { g_i3->b_subvertex.v->subvertex = g_i3; }
576 if (g_i2_b_subvertex_dim == 0 && g_i2 != g_i3) { b_edge->subvertex = g_i3; }
577 if (g_i3->b_subvertex.dim == 0 && g_i3->b_subvertex.v == b_edge->vertex) {
578 if (g_i2 != g_i3) {
579 b_edge->pair->subvertex = g_i2;
580 } else if (g_i1) {
581 b_edge->pair->subvertex = g_i1;
582 }
583 b_edge->vertex->subvertex = g_i3;
584 break;
585 }
586 if (g_i3 != g_i2) {
587 g_i1 = g_i2;
588 g_i2 = g_i3;
589 }
590 }
591 if (g_i3) {
592 b_vertex_iterator_on_edge iter(b_edge->vertex);
593 for (b_edge_type* e = iter.next(); e != 0; e = iter.next()) {
594 edge_to_process.push(e);
595 }
596 }
597 }
598 }
599
600 // step 1 twice: all the subvertex on b and g edge must point on the lower edge
601 for (typename list_subvertex_t::iterator sv = list_subvertex.begin();
602 sv != list_subvertex.end(); ++sv) {
603 if (sv->b_subvertex.dim == 1) {
604 b_edge_type* e = sv->b_subvertex.e;
605 if (e > e->pair) {
606 sv->b_subvertex.e = e->pair;
607 sv->b_subvertex.icoor[0] = 1 - sv->b_subvertex.icoor[0];
608 }
609 }
610 if (sv->g_subvertex.dim == 1) {
611 g_edge_type* e = sv->g_subvertex.e;
612 if (e > e->pair) {
613 sv->g_subvertex.e = e->pair;
614 sv->g_subvertex.icoor[0] = 1 - sv->g_subvertex.icoor[0];
615 }
616 }
617 }
618
619#ifdef B2OSM_DEBUG
620 std::cout << "B VERTEX" << std::endl;
621 for (size_t i = 0; i != b_mesh->list_vertex.size(); ++i) {
622 b_vertex_type& v = b_mesh->list_vertex[i];
623 std::cout << "id = " << v.id << ", sv = " << get_subvertex_id(v.subvertex) << std::endl;
624 }
625
626 std::cout << "B EDGE" << std::endl;
627 for (size_t i = 0; i != b_mesh->list_edge.size(); ++i) {
628 b_edge_type& e = b_mesh->list_edge[i];
629 SubVertex* sv = static_cast<SubVertex*>(e.subvertex);
630 std::cout << "id = " << e.pair->vertex->id << ", " << e.vertex->id
631 << ", sv = " << get_subvertex_id(sv);
632 if (sv) {
633 std::cout << ", sv dim = " << sv->b_subvertex.dim
634 << ", sv icoor = " << sv->b_subvertex.icoor[0];
635 }
636 std::cout << std::endl;
637 }
638
639 std::cout << "G VERTEX" << std::endl;
640 for (size_t i = 0; i != g_mesh->list_vertex.size(); ++i) {
641 g_vertex_type& v = g_mesh->list_vertex[i];
642 std::cout << "id = " << v.id << ", sv = " << get_subvertex_id(v.subvertex) << std::endl;
643 }
644
645 std::cout << "G EDGE" << std::endl;
646 for (size_t i = 0; i != g_mesh->list_edge.size(); ++i) {
647 g_edge_type& e = g_mesh->list_edge[i];
648 std::cout << "id = " << e.pair->vertex->id << ", " << e.vertex->id
649 << ", sv = " << get_subvertex_id(e.subvertex) << std::endl;
650 }
651
652 std::cout << "SUBVERTEX" << std::endl;
653 for (size_t i = 0; i != list_subvertex.size(); ++i) { print_subvertex(list_subvertex[i]); }
654 std::cout << std::endl;
655#endif
656
657 // step two: sorting subvertices in green edge
658 {
659 typedef std::list<std::pair<g_edge_type*, SubVertex*> > list_sv_t;
660 list_sv_t list_sv;
661 // we iterate on the blue face
662 for (typename B_DCEL_MESH::face_iterator bf = b_mesh->face_begin();
663 bf != b_mesh->face_end(); ++bf) {
664 list_sv.clear();
665 b_edge_type* b_edge = bf->edge;
666 b_edge_type* b_edge_start = b_edge;
667 do {
668 SubVertex* sv = static_cast<SubVertex*>(b_edge->subvertex);
669 for (;;) {
670 if (sv->g_subvertex.dim == 1) {
671 g_edge_type* e = sv->g_subvertex.e;
672 if (e > e->pair) { e = e->pair; }
673 list_sv.push_back(std::pair<g_edge_type*, SubVertex*>(e, sv));
674 } else if (sv->g_subvertex.dim == 0) {
675 g_vertex_iterator_on_edge iter(sv->g_subvertex.v);
676 for (g_edge_type* e = iter.next(); e != 0; e = iter.next()) {
677 list_sv.push_back(std::pair<g_edge_type*, SubVertex*>(
678 e->pair < e ? e->pair : e, sv));
679 }
680 }
681 if (sv->b_subvertex.dim == 0) { break; }
682 if (b_edge < b_edge->pair) {
683 sv = sv->b_next_subvertex;
684 } else {
685 sv = sv->b_previous_subvertex;
686 }
687 }
688 b_edge = b_edge->next;
689 } while (b_edge != b_edge_start);
690
691 if (list_sv.size() < 2) { continue; }
692
693 list_sv.sort(CompSV());
694
695#ifdef B2OSM_DEBUG
696 {
697 std::cout << "for face " << bf->id << std::endl;
698 for (typename list_sv_t::iterator i = list_sv.begin(); i != list_sv.end();
699 ++i) {
700 std::cout << i->first << ", ";
701 print_subvertex(*i->second);
702 }
703 }
704#endif
705
706 // Remove extra subvertex
707 {
708 typename list_sv_t::iterator i_in = list_sv.begin();
709
710 while (i_in != list_sv.end()) {
711 typename list_sv_t::iterator i_in1 = i_in;
712 size_t s = 1;
713 ++i_in1;
714 while (i_in1 != list_sv.end() && i_in1->first == i_in->first) {
715 ++i_in1;
716 ++s;
717 }
718 if (s > 2 && bf != b_mesh->face_end()) {
719 typedef std::deque<std::pair<double, typename list_sv_t::iterator> >
720 list_diff_t;
721 list_diff_t list_diff;
722 for (typename list_sv_t::iterator i = i_in; i != i_in1; ++i) {
723 if (i->second->g_subvertex.dim == 0) {
724 if (i->second->g_subvertex.v == i->first->vertex) {
725 list_diff.push_back(typename list_diff_t::value_type(1, i));
726 } else {
727 list_diff.push_back(typename list_diff_t::value_type(0, i));
728 }
729 } else {
730 list_diff.push_back(typename list_diff_t::value_type(
731 i->second->g_subvertex.icoor[0], i));
732 }
733 }
734 for (size_t i = 0; i != list_diff.size() - 1; ++i) {
735 list_diff[i].first = list_diff[i + 1].first - list_diff[i].first;
736 }
737 list_diff.pop_back();
738 std::sort(list_diff.begin(), list_diff.end(), CompareFirstOfPair());
739
740 while (list_diff.size() > 1) {
741 typename list_sv_t::iterator a = list_diff.front().second;
742 typename list_sv_t::iterator b = a;
743 ++b;
744 cluster_subvertex_on_dcel(a->second, b->second);
745 while (i_in1 != list_sv.end() && i_in1->second == b->second) {
746 ++i_in1;
747 }
748 {
749 SubVertex* v = b->second;
750 for (typename list_sv_t::iterator i = i_in;
751 i != list_sv.end();) {
752 if (i->second == v) {
753 i = list_sv.erase(i);
754 } else {
755 ++i;
756 }
757 }
758 }
759 list_diff.pop_front();
760 }
761 }
762 i_in = i_in1;
763 }
764 }
765
766#ifdef B2OSM_DEBUG
767 {
768 std::cout << "after remove extra subvertex on face " << bf->id << std::endl;
769 for (typename list_sv_t::iterator i = list_sv.begin(); i != list_sv.end();
770 ++i) {
771 std::cout << i->first << ", ";
772 print_subvertex(*i->second);
773 }
774 }
775#endif
776
777 typename list_sv_t::iterator i_end = list_sv.end();
778 --i_end;
779 for (typename list_sv_t::iterator i = list_sv.begin(); i != i_end; ++i) {
780 typename list_sv_t::iterator j = i;
781 ++j;
782 if (j->first == i->first) {
783 SubVertex* sv1 = j->second;
784 SubVertex* sv2 = i->second;
785 if (sv1->g_subvertex.dim == 1 && sv2->g_subvertex.dim == 1) {
786 if (sv1->g_subvertex.icoor[0] > sv2->g_subvertex.icoor[0]) {
787 std::swap(sv1, sv2);
788 }
789 } else if (sv1->g_subvertex.dim == 1 || sv2->g_subvertex.dim == 1) {
790 if (sv1->g_subvertex.dim == 1) { std::swap(sv1, sv2); }
791 if (sv2->g_subvertex.e->vertex == sv1->g_subvertex.v) {
792 std::swap(sv1, sv2);
793 }
794 }
795 sv1->g_next_subvertex = sv2;
796 sv2->g_previous_subvertex = sv1;
797 }
798 }
799 }
800 }
801
802#ifdef B2OSM_DEBUG
803 std::cout << "B VERTEX" << std::endl;
804 for (size_t i = 0; i != b_mesh->list_vertex.size(); ++i) {
805 b_vertex_type& v = b_mesh->list_vertex[i];
806 std::cout << "id = " << v.id << ", sv = " << get_subvertex_id(v.subvertex) << std::endl;
807 }
808
809 std::cout << "B EDGE" << std::endl;
810 for (size_t i = 0; i != b_mesh->list_edge.size(); ++i) {
811 b_edge_type& e = b_mesh->list_edge[i];
812 SubVertex* sv = static_cast<SubVertex*>(e.subvertex);
813 std::cout << "id = " << e.pair->vertex->id << ", " << e.vertex->id
814 << ", sv = " << get_subvertex_id(sv);
815 if (sv) {
816 std::cout << ", sv dim = " << sv->b_subvertex.dim
817 << ", sv icoor = " << sv->b_subvertex.icoor[0];
818 }
819 std::cout << std::endl;
820 }
821
822 std::cout << "G VERTEX" << std::endl;
823 for (size_t i = 0; i != g_mesh->list_vertex.size(); ++i) {
824 g_vertex_type& v = g_mesh->list_vertex[i];
825 std::cout << "id = " << v.id << ", sv = " << get_subvertex_id(v.subvertex) << std::endl;
826 }
827
828 std::cout << "G EDGE" << std::endl;
829 for (size_t i = 0; i != g_mesh->list_edge.size(); ++i) {
830 g_edge_type& e = g_mesh->list_edge[i];
831 std::cout << "id = " << e.pair->vertex->id << ", " << e.vertex->id
832 << ", sv = " << get_subvertex_id(e.subvertex) << std::endl;
833 }
834
835 std::cout << "SUBVERTEX" << std::endl;
836 for (size_t i = 0; i != list_subvertex.size(); ++i) { print_subvertex(list_subvertex[i]); }
837 std::cout << std::endl;
838#endif
839
840 // step three: projecting green vertices
841 for (typename B_DCEL_MESH::face_iterator bf = b_mesh->face_begin();
842 bf != b_mesh->face_end(); ++bf) {
843 std::vector<void*> list_b_cell;
844 list_b_cell.push_back(&*bf);
845 b_edge_type* be_start;
846 be_start = bf->edge;
847
848 b_edge_type* be = be_start;
849 do {
850 list_b_cell.push_back(be);
851 list_b_cell.push_back(be->vertex);
852 be = be->next;
853 } while (be != be_start);
854 std::sort(list_b_cell.begin(), list_b_cell.end());
855
856 be = be_start;
857 do {
858 SubVertex* sv = static_cast<SubVertex*>(be->subvertex);
859 for (;;) {
860 if (sv->g_subvertex.dim == 0) {
861 g_vertex_iterator_on_edge iter(sv->g_subvertex.v);
862 for (g_edge_type* ge = iter.next(); ge != 0; ge = iter.next()) {
863 if (!ge->subvertex) {
864 try_g_vertice_on_blue_face(&*bf, list_b_cell, ge);
865 }
866 }
867 } else if (
868 sv->g_subvertex.dim == 1
869 && !((sv->g_previous_subvertex && sv->g_next_subvertex)
870 || (sv->g_previous_subvertex
871 && std::binary_search(
872 list_b_cell.begin(), list_b_cell.end(),
873 sv->g_previous_subvertex->g_subvertex.v))
874 || (sv->g_next_subvertex
875 && std::binary_search(
876 list_b_cell.begin(), list_b_cell.end(),
877 sv->g_next_subvertex->g_subvertex.v)))) {
878 bool res = false;
879 if (!sv->g_next_subvertex) {
880 res = try_g_vertice_on_blue_face(&*bf, list_b_cell, sv->g_subvertex.e);
881 }
882 if (!res && sv->g_previous_subvertex) {
883 try_g_vertice_on_blue_face(&*bf, list_b_cell, sv->g_subvertex.e->pair);
884 }
885 }
886
887 if (sv->b_subvertex.dim == 1) {
888 if (sv->b_subvertex.e == be) {
889 sv = sv->b_next_subvertex;
890 } else {
891 assert(sv->b_subvertex.e->pair == be);
892 sv = sv->b_previous_subvertex;
893 }
894 } else {
895 break;
896 }
897 }
898 be = be->next;
899 } while (be != be_start);
900 }
901
902#ifdef B2OSM_DEBUG
903 std::cout << "B VERTEX" << std::endl;
904 for (size_t i = 0; i != b_mesh->list_vertex.size(); ++i) {
905 b_vertex_type& v = b_mesh->list_vertex[i];
906 std::cout << "id = " << v.id << ", sv = " << get_subvertex_id(v.subvertex) << std::endl;
907 }
908
909 std::cout << "B EDGE" << std::endl;
910 for (size_t i = 0; i != b_mesh->list_edge.size(); ++i) {
911 b_edge_type& e = b_mesh->list_edge[i];
912 SubVertex* sv = static_cast<SubVertex*>(e.subvertex);
913 std::cout << "id = " << e.pair->vertex->id << ", " << e.vertex->id
914 << ", sv = " << get_subvertex_id(sv);
915 if (sv) {
916 std::cout << ", sv dim = " << sv->b_subvertex.dim
917 << ", sv icoor = " << sv->b_subvertex.icoor[0];
918 }
919 std::cout << std::endl;
920 }
921
922 std::cout << "G VERTEX" << std::endl;
923 for (size_t i = 0; i != g_mesh->list_vertex.size(); ++i) {
924 g_vertex_type& v = g_mesh->list_vertex[i];
925 std::cout << "id = " << v.id << ", sv = " << get_subvertex_id(v.subvertex) << std::endl;
926 }
927
928 std::cout << "G EDGE" << std::endl;
929 for (size_t i = 0; i != g_mesh->list_edge.size(); ++i) {
930 g_edge_type& e = g_mesh->list_edge[i];
931 std::cout << "id = " << e.pair->vertex->id << ", " << e.vertex->id
932 << ", sv = " << get_subvertex_id(e.subvertex) << std::endl;
933 }
934
935 std::cout << "SUBVERTEX" << std::endl;
936 for (size_t i = 0; i != list_subvertex.size(); ++i) { print_subvertex(list_subvertex[i]); }
937 std::cout << std::endl;
938#endif
939
940 for (typename B_DCEL_MESH::vertex_iterator v = b_mesh->vertex_begin();
941 v != b_mesh->vertex_end(); ++v) {
942 if (v->subvertex && static_cast<SubVertex*>(v->subvertex)->g_subvertex.dim == 0) {}
943 }
944
945 // stage 3 twice:
946 // - set green vertices that have no subvertex has subvertex that do not
947 // project on the blue mesh,
948 // - set pointer to the subvertex on the green edges,
949 // - connect the subvertex on green edge with the subvertex on green vertex.
950 for (typename G_DCEL_MESH::vertex_iterator v = g_mesh->vertex_begin();
951 v != g_mesh->vertex_end(); ++v) {
952 if (v->subvertex == 0) {
953 SubVertex& sv = new_sub_vertex();
954 v->subvertex = &sv;
955 sv.b_subvertex.dim = 2;
956 sv.b_subvertex.f = 0;
957 std::fill_n(sv.b_subvertex.icoor, 2, 0.0);
958 sv.g_subvertex.dim = 0;
959 sv.g_subvertex.v = &*v;
960 }
961 }
962 for (typename list_subvertex_t::iterator sv = list_subvertex.begin();
963 sv != list_subvertex.end(); ++sv) {
964 if (sv->g_subvertex.dim == 1) {
965 if (sv->g_next_subvertex == 0) {
966 sv->g_subvertex.e->pair->subvertex = &*sv;
967 sv->g_next_subvertex =
968 static_cast<SubVertex*>(sv->g_subvertex.e->vertex->subvertex);
969 sv->g_next_subvertex->g_previous_subvertex = &*sv;
970 } else if (sv->g_next_subvertex->g_subvertex.dim == 0) {
971 sv->g_subvertex.e->pair->subvertex = &*sv;
972 }
973 if (sv->g_previous_subvertex == 0) {
974 sv->g_subvertex.e->subvertex = &*sv;
975 sv->g_previous_subvertex =
976 static_cast<SubVertex*>(sv->g_subvertex.e->pair->vertex->subvertex);
977 sv->g_previous_subvertex->g_next_subvertex = &*sv;
978 } else if (sv->g_previous_subvertex->g_subvertex.dim == 0) {
979 sv->g_subvertex.e->subvertex = &*sv;
980 }
981 }
982 }
983 for (typename G_DCEL_MESH::edge_iterator e = g_mesh->edge_begin(); e != g_mesh->edge_end();
984 ++e) {
985 if (e->subvertex == 0) { e->subvertex = e->vertex->subvertex; }
986 }
987
988#ifdef B2OSM_DEBUG
989 std::cout << "B VERTEX" << std::endl;
990 for (size_t i = 0; i != b_mesh->list_vertex.size(); ++i) {
991 b_vertex_type& v = b_mesh->list_vertex[i];
992 std::cout << "id = " << v.id << ", sv = " << get_subvertex_id(v.subvertex) << std::endl;
993 }
994
995 std::cout << "B EDGE" << std::endl;
996 for (size_t i = 0; i != b_mesh->list_edge.size(); ++i) {
997 b_edge_type& e = b_mesh->list_edge[i];
998 SubVertex* sv = static_cast<SubVertex*>(e.subvertex);
999 std::cout << "id = " << e.pair->vertex->id << ", " << e.vertex->id << ", ptr = " << &e
1000 << ", sv = " << get_subvertex_id(sv);
1001 if (sv) {
1002 std::cout << ", sv dim = " << sv->b_subvertex.dim
1003 << ", sv icoor = " << sv->b_subvertex.icoor[0];
1004 }
1005 std::cout << std::endl;
1006 }
1007
1008 std::cout << "G VERTEX" << std::endl;
1009 for (size_t i = 0; i != g_mesh->list_vertex.size(); ++i) {
1010 g_vertex_type& v = g_mesh->list_vertex[i];
1011 std::cout << "id = " << v.id << ", sv = " << get_subvertex_id(v.subvertex) << std::endl;
1012 }
1013
1014 std::cout << "G EDGE" << std::endl;
1015 for (size_t i = 0; i != g_mesh->list_edge.size(); ++i) {
1016 g_edge_type& e = g_mesh->list_edge[i];
1017 std::cout << "id = " << e.pair->vertex->id << ", " << e.vertex->id << ", ptr = " << &e
1018 << ", sv = " << get_subvertex_id(e.subvertex) << std::endl;
1019 }
1020
1021 std::cout << "SUBVERTEX" << std::endl;
1022 for (size_t i = 0; i != list_subvertex.size(); ++i) { print_subvertex(list_subvertex[i]); }
1023 std::cout << std::endl;
1024
1025 std::cout << "branch 3" << std::endl;
1026 std::cout << "nodes" << std::endl;
1027 for (size_t i = 0; i != list_subvertex.size(); ++i) {
1028 std::cout << i + 1 << " ";
1029 print_subvertex_bcoor(list_subvertex[i]);
1030 std::cout << std::endl;
1031 }
1032 std::cout << "end\nelements\n type PMASS3\n mass_coefficients_3 0.1 0.1 0.1\n1 "
1033 "1\nend\nendbranch"
1034 << std::endl;
1035
1036 std::cout << "branch 4" << std::endl;
1037 std::cout << "nodes" << std::endl;
1038 for (size_t i = 0; i != list_subvertex.size(); ++i) {
1039 std::cout << i + 1 << " ";
1040 print_subvertex_gcoor(list_subvertex[i]);
1041 std::cout << std::endl;
1042 }
1043 std::cout << "end\nelements\n type PMASS3\n mass_coefficients_3 0.1 0.1 0.1\n1 "
1044 "1\nend\nendbranch"
1045 << std::endl;
1046#endif
1047 }
1048
1054 class OverlayingMeshIterator {
1055 public:
1056 typedef ComputeOverlayingSurfaceMeshes<G_DCEL_MESH, B_DCEL_MESH> parent_type;
1057 OverlayingMeshIterator(
1058 const ComputeOverlayingSurfaceMeshes& parent_, bool non_overlapping_triangle_)
1059 : parent(parent_),
1060 b_face(parent_.b_mesh->face_begin()),
1061 b_face_end(parent_.b_mesh->face_end()),
1062 non_overlapping_triangle(non_overlapping_triangle_) {}
1063
1064 template <typename TRIANGLE>
1065 bool next(TRIANGLE& triangle) {
1066 for (;;) {
1067 if (b_face == b_face_end && !non_overlapping_triangle) { return false; }
1068
1069 list_in_edge_of_b_face_t list_sv;
1070 if (stack.empty()) {
1071 // we initialise the list of sub-vertices in the edge of the
1072 // current blue face
1073 // list_sv.reserve(3);
1074 b_edge_type* b_edge;
1075 if (b_face == b_face_end) {
1076 b_edge = parent.b_border_edge;
1077 } else {
1078 b_edge = b_face->edge;
1079 }
1080 b_edge_type* b_edge_end = b_edge;
1081 do {
1082 typename parent_type::SubVertex* v =
1083 static_cast<typename parent_type::SubVertex*>(b_edge->subvertex);
1084 while (v != 0) {
1085 list_sv.push_back(v);
1086 if (v->b_subvertex.dim == 0) { break; }
1087 if (b_face == b_face_end) {
1088 if (b_edge < b_edge->pair) {
1089 v = v->b_previous_subvertex;
1090 } else {
1091 v = v->b_next_subvertex;
1092 }
1093 } else if (b_edge < b_edge->pair) {
1094 v = v->b_next_subvertex;
1095 } else {
1096 v = v->b_previous_subvertex;
1097 }
1098 }
1099 b_edge = b_edge->next;
1100 } while (b_edge != b_edge_end);
1101 } else {
1102 list_sv.swap(stack.top());
1103 stack.pop();
1104 }
1105
1106#ifdef B2OSM_DEBUG
1107 std::cout << "SUBV start" << std::endl;
1108 for (typename list_in_edge_of_b_face_t::iterator i = list_sv.begin();
1109 i != list_sv.end(); ++i) {
1110 parent.print_subvertex(**i);
1111 }
1112 std::cout << std::endl;
1113#endif
1114
1115 triangle.b_face = b_face == b_face_end ? 0 : &*b_face;
1116 typename list_in_edge_of_b_face_t::iterator sv1 = list_sv.begin();
1117 typename list_in_edge_of_b_face_t::iterator sv2 = sv1;
1118 ++sv2;
1119 bool list_sv_non_overlapping_polygon = false;
1120 try {
1121 triangle.g_face = get_face((*sv1)->g_subvertex, (*sv2)->g_subvertex);
1122 } catch (std::exception& e) { list_sv_non_overlapping_polygon = true; }
1123 if (!list_sv_non_overlapping_polygon && triangle.g_face == 0) {
1124 typename parent_type::SubVertex* sv1_value = *sv1;
1125 for (;;) {
1126 list_sv.splice(list_sv.end(), list_sv, sv1);
1127 sv1 = list_sv.begin();
1128 if (sv1_value == *sv1) { break; }
1129 sv2 = sv1;
1130 ++sv2;
1131 triangle.g_face = get_face((*sv1)->g_subvertex, (*sv2)->g_subvertex);
1132 if (triangle.g_face != 0) { break; }
1133 }
1134 if (sv1_value == *sv1) {
1135 // The list_sv define a non overlapping polygon
1136 if (!non_overlapping_triangle) {
1137 list_sv_non_overlapping_polygon = true;
1138 } else {
1139 UnimplementedError() << THROW;
1140 }
1141 }
1142 }
1143
1144 if (!list_sv_non_overlapping_polygon) {
1145 SubVertex* sv3 = 0;
1146 try {
1147 sv3 = get_third_subvertex(*sv1, *sv2, triangle.g_face, triangle.b_face);
1148 } catch (std::exception& e) { ; }
1149 if (sv3) {
1150#ifdef B2OSM_DEBUG
1151 std::cout << "1 :";
1152 parent.print_subvertex(**sv1);
1153 std::cout << "2 :";
1154 parent.print_subvertex(**sv2);
1155 std::cout << "3 :";
1156 parent.print_subvertex(*sv3);
1157 std::cout << std::endl;
1158#endif
1159
1160 set_coor_in_triangle(triangle, *sv1, 0);
1161 set_coor_in_triangle(triangle, *sv2, 1);
1162 set_coor_in_triangle(triangle, sv3, 2);
1163
1164 // check if the sv3 subvertex is in the list_sv
1165 typename list_in_edge_of_b_face_t::iterator sv3_iter =
1166 std::find(sv2, list_sv.end(), sv3);
1167 if (sv3_iter == list_sv.end()) {
1168 // insert sv3 between sv1 and sv2 and push it in the stack
1169 list_sv.insert(sv2, sv3);
1170 stack.push(list_in_edge_of_b_face_t());
1171 list_sv.swap(stack.top());
1172 } else {
1173 // the sv2->sv3 or sv3->sv1 edges split the list_sv in two list
1174 typename list_in_edge_of_b_face_t::iterator tmp = sv3_iter;
1175 if (tmp == list_sv.begin()) { tmp = list_sv.end(); }
1176 --tmp;
1177 if (tmp != sv2) {
1178 stack.push(list_in_edge_of_b_face_t());
1179 stack.top().splice(stack.top().end(), list_sv, sv2, sv3_iter);
1180 stack.top().push_back(sv3);
1181
1182#ifdef B2OSM_DEBUG
1183 std::cout << "SUBV push" << std::endl;
1184 for (typename list_in_edge_of_b_face_t::iterator i =
1185 list_sv.begin();
1186 i != list_sv.end(); ++i) {
1187 parent.print_subvertex(**i);
1188 }
1189 std::cout << std::endl;
1190#endif
1191 } else {
1192 list_sv.erase(sv2);
1193 }
1194 tmp = sv3_iter;
1195 ++tmp;
1196 if (tmp == list_sv.end()) { tmp = list_sv.begin(); }
1197 if (tmp != sv1) {
1198 stack.push(list_in_edge_of_b_face_t());
1199 list_sv.swap(stack.top());
1200 }
1201 }
1202 }
1203 }
1204
1205 // if the stack of list of subvertex to process is empty, we pass to
1206 // the next blue face.
1207 if (stack.empty()) {
1208 if (b_face == b_face_end) {
1209 assert(non_overlapping_triangle == true);
1210 non_overlapping_triangle = false;
1211 if (triangle.b_face != 0 || triangle.g_face != 0) { return true; }
1212 } else {
1213 ++b_face;
1214 }
1215 }
1216 if ((non_overlapping_triangle && (triangle.b_face != 0 || triangle.g_face != 0))
1217 || (triangle.b_face != 0 && triangle.g_face != 0)) {
1218 return true;
1219 }
1220 }
1221 }
1222
1223 private:
1224 template <typename TRIANGLE>
1225 static void set_coor_in_triangle(
1226 TRIANGLE& triangle, typename parent_type::SubVertex* sv, int coor_id) {
1227 if (triangle.g_face != 0) {
1228 get_internal_coor(
1229 *triangle.g_face, sv->g_subvertex, triangle.internal_coor_g[coor_id]);
1230 }
1231 if (triangle.b_face != 0) {
1232 get_internal_coor(
1233 *triangle.b_face, sv->b_subvertex, triangle.internal_coor_b[coor_id]);
1234 }
1235 };
1236
1237 template <typename FACE_TYPE, typename XSUBVERTEX>
1238 static void get_internal_coor(
1239 const FACE_TYPE& face, const XSUBVERTEX& sv, double icoor[2]) {
1240 switch (sv.dim) {
1241 case 0: {
1242 typename FACE_TYPE::edge_type* e = face.edge;
1243 if (e->vertex == sv.v) {
1244 icoor[0] = icoor[1] = 0;
1245 return;
1246 }
1247 e = e->next;
1248 if (e->vertex == sv.v) {
1249 icoor[0] = 1;
1250 icoor[1] = 0;
1251 return;
1252 }
1253 assert(e->next->vertex == sv.v);
1254 icoor[0] = 0;
1255 icoor[1] = 1;
1256 return;
1257 }
1258 case 1: {
1259 typename FACE_TYPE::edge_type* e = face.edge;
1260 typename FACE_TYPE::edge_type* sve = sv.e < sv.e->pair ? sv.e : sv.e->pair;
1261 const double c = sv.icoor[0];
1262 if (e == sve) {
1263 icoor[0] = 0;
1264 icoor[1] = 1 - c;
1265 return;
1266 }
1267 if (e == sve->pair) {
1268 icoor[0] = 0;
1269 icoor[1] = c;
1270 return;
1271 }
1272 e = e->next;
1273 if (e == sve) {
1274 icoor[0] = c;
1275 icoor[1] = 0;
1276 return;
1277 }
1278 if (e == sve->pair) {
1279 icoor[0] = 1 - c;
1280 icoor[1] = 0;
1281 return;
1282 }
1283 e = e->next;
1284 if (e == sve) {
1285 icoor[0] = 1 - c;
1286 icoor[1] = c;
1287 return;
1288 }
1289 assert(e == sve->pair);
1290 icoor[0] = c;
1291 icoor[1] = 1 - c;
1292 return;
1293 }
1294 case 2:
1295 std::copy(sv.icoor, sv.icoor + 2, icoor);
1296 return;
1297 }
1298 }
1299
1300 template <typename SUB_VERTEX>
1301 static typename SUB_VERTEX::dcel_mesh_type::face_type* get_face(
1302 SUB_VERTEX& sv1, SUB_VERTEX& sv2) {
1303 typedef typename SUB_VERTEX::dcel_mesh_type::face_type face_type;
1304 if (sv1.dim == 2) { return sv1.f; }
1305 if (sv2.dim == 2) { return sv2.f; }
1306 if (sv1.dim == 1 && sv2.dim == 1) {
1307 if (sv1.e == sv2.e) {
1308 if (sv1.icoor[0] < sv2.icoor[0]) {
1309 return sv1.e->face;
1310 } else {
1311 return sv1.e->pair->face;
1312 }
1313 }
1314 face_type* f11 = sv1.e->face;
1315 face_type* f21 = sv2.e->face;
1316 if (f11 == f21) { return f11; }
1317 if (sv2.e->pair->face == 0) { return f21; }
1318 face_type* f12 = sv1.e->pair->face;
1319 if (f12 == f21) { return f12; }
1320 assert(sv2.e->pair->face != 0);
1321 return sv2.e->pair->face;
1322 }
1323 if (sv1.dim == 0 && sv2.dim == 0) {
1324 typename SUB_VERTEX::dcel_mesh_type::edge_type* e = sv1.v->edge;
1325 do {
1326 if (e->vertex == sv2.v) { return e->face; }
1327 e = e->pair->next;
1328 } while (e != sv1.v->edge);
1329 e = sv2.v->edge;
1330 do {
1331 if (e->vertex == sv1.v) { return e->pair->face; }
1332 e = e->pair->next;
1333 } while (e != sv2.v->edge);
1334
1335 {
1336 std::vector<void*> sv1_face;
1337 sv1_face.reserve(10);
1338 e = sv1.v->edge;
1339 do {
1340 sv1_face.push_back(e->face);
1341 e = e->pair->next;
1342 } while (e != sv1.v->edge);
1343 std::sort(sv1_face.begin(), sv1_face.end());
1344 e = sv2.v->edge;
1345 do {
1346 if (std::find(sv1_face.begin(), sv1_face.end(), e->face)
1347 != sv1_face.end()) {
1348 return e->face;
1349 }
1350 e = e->pair->next;
1351 } while (e != sv2.v->edge);
1352 }
1353 assert(false);
1354 }
1355 if (sv1.dim == 0) {
1356 typename SUB_VERTEX::dcel_mesh_type::edge_type* e = sv2.e;
1357 do {
1358 if (e->vertex == sv1.v) {
1359 if (e == sv2.e) { e = e->pair; }
1360 return e->face;
1361 }
1362 e = e->next;
1363 } while (e != sv2.e);
1364 e = sv2.e->pair;
1365 do {
1366 if (e->vertex == sv1.v) {
1367 if (e == sv2.e->pair) { e = e->pair; }
1368 return e->face;
1369 }
1370 e = e->next;
1371 } while (e != sv2.e->pair);
1372 } else {
1373 typename SUB_VERTEX::dcel_mesh_type::edge_type* e = sv1.e;
1374 do {
1375 if (e->vertex == sv2.v) {
1376 if (e->next == sv1.e) { e = e->next->pair; }
1377 return e->face;
1378 }
1379 e = e->next;
1380 } while (e != sv1.e);
1381 e = sv1.e->pair;
1382 do {
1383 if (e->vertex == sv2.v) {
1384 if (e->next == sv1.e->pair) { e = e->next->pair; }
1385 return e->face;
1386 }
1387 e = e->next;
1388 } while (e != sv1.e->pair);
1389 }
1390#ifndef AVOID_PROBLEM_IN_TRIANGLE_ITERATOR
1391 assert(false);
1392#endif
1393 throw std::exception();
1394 return 0;
1395 }
1396
1397 static typename parent_type::SubVertex* get_third_subvertex(
1398 typename parent_type::SubVertex* sv1, typename parent_type::SubVertex* sv2,
1399 typename parent_type::g_face_type* g_face,
1400 typename parent_type::b_face_type* b_face) {
1401 SubVertex* sv3 = 0;
1402 if (sv2->b_subvertex.dim == 1) {
1403 sv3 = sv2->b_next_subvertex;
1404 if (sv3 != sv1 && get_face(sv2->g_subvertex, sv3->g_subvertex) == g_face
1405 && get_face(sv2->b_subvertex, sv3->b_subvertex) == b_face) {
1406 return sv3;
1407 }
1408 sv3 = sv2->b_previous_subvertex;
1409 if (sv3 != sv1 && get_face(sv2->g_subvertex, sv3->g_subvertex) == g_face
1410 && get_face(sv2->b_subvertex, sv3->b_subvertex) == b_face) {
1411 return sv3;
1412 }
1413 }
1414 if (sv2->g_subvertex.dim == 1) {
1415 sv3 = sv2->g_next_subvertex;
1416 if (sv3 != sv1 && get_face(sv2->g_subvertex, sv3->g_subvertex) == g_face
1417 && get_face(sv2->b_subvertex, sv3->b_subvertex) == b_face) {
1418 return sv3;
1419 }
1420 sv3 = sv2->g_previous_subvertex;
1421 if (sv3 != sv1 && get_face(sv2->g_subvertex, sv3->g_subvertex) == g_face
1422 && get_face(sv2->b_subvertex, sv3->b_subvertex) == b_face) {
1423 return sv3;
1424 }
1425 }
1426
1427 if (sv2->b_subvertex.dim == 0) {
1428 b_vertex_iterator_on_edge iter(sv2->b_subvertex.v);
1429 for (b_edge_type* e = iter.next(); e != 0; e = iter.next()) {
1430 sv3 = static_cast<SubVertex*>(e->subvertex);
1431 if (sv3 != sv1 && get_face(sv2->g_subvertex, sv3->g_subvertex) == g_face
1432 && get_face(sv2->b_subvertex, sv3->b_subvertex) == b_face) {
1433 return sv3;
1434 }
1435 }
1436 }
1437
1438 if (sv2->g_subvertex.dim == 0) {
1439 g_vertex_iterator_on_edge iter(sv2->g_subvertex.v);
1440 for (g_edge_type* e = iter.next(); e != 0; e = iter.next()) {
1441 sv3 = static_cast<SubVertex*>(e->subvertex);
1442 if (sv3 != sv1 && get_face(sv2->g_subvertex, sv3->g_subvertex) == g_face
1443 && get_face(sv2->b_subvertex, sv3->b_subvertex) == b_face) {
1444 return sv3;
1445 }
1446 }
1447 }
1448#ifndef AVOID_PROBLEM_IN_TRIANGLE_ITERATOR
1449 assert(false);
1450#endif
1451 return 0;
1452 }
1453
1454 const ComputeOverlayingSurfaceMeshes& parent;
1455 typename B_DCEL_MESH::face_iterator b_face;
1456 typename B_DCEL_MESH::face_iterator b_face_end;
1457 typedef std::list<typename parent_type::SubVertex*> list_in_edge_of_b_face_t;
1458 typedef std::stack<list_in_edge_of_b_face_t> stack_t;
1459 stack_t stack;
1460 bool non_overlapping_triangle;
1461 };
1462
1463 OverlayingMeshIterator get_overlaying_mesh_iterator(bool non_overlapping_triangle = false) {
1464 return OverlayingMeshIterator(*this, non_overlapping_triangle);
1465 }
1466
1467private:
1468 friend class OverlayingMeshIterator;
1469
1470 class b_vertex_iterator_on_edge {
1471 public:
1472 b_vertex_iterator_on_edge(b_vertex_type* v_) : v(v_), e(v_->edge) {}
1473
1474 b_edge_type* next() {
1475 b_edge_type* res = e;
1476 if (e != 0) {
1477 e = e->pair->next;
1478 if (e == v->edge) { e = 0; }
1479 }
1480 return res;
1481 }
1482
1483 private:
1484 b_vertex_type* v;
1485 b_edge_type* e;
1486 };
1487
1488 class g_vertex_iterator_on_edge {
1489 public:
1490 g_vertex_iterator_on_edge(g_vertex_type* v_) : v(v_), e(v_->edge) {}
1491
1492 g_edge_type* next() {
1493 g_edge_type* res = e;
1494 if (e != 0) {
1495 e = e->pair->next;
1496 if (e == v->edge) { e = 0; }
1497 }
1498 return res;
1499 }
1500
1501 private:
1502 g_vertex_type* v;
1503 g_edge_type* e;
1504 };
1505
1506 struct GCell {
1507 int dim; // 0 is vertex, 1 is edge, 2 is face, -1 is none
1508 union {
1509 g_vertex_type* v;
1510 g_edge_type* e;
1511 g_face_type* f;
1512 };
1513 bool operator<(const GCell& c) const { return v < c.v; }
1514 };
1515
1516 struct BCell {
1517 int dim; // 0 is vertex, 1 is edge, 2 is face, -1 is none
1518 union {
1519 b_vertex_type* v;
1520 b_edge_type* e;
1521 b_face_type* f;
1522 };
1523 bool operator<(const BCell& c) const { return v < c.v; }
1524 };
1525
1526 struct GSubVertex : public GCell {
1527 typedef G_DCEL_MESH dcel_mesh_type;
1528 double icoor[2];
1529 };
1530
1531 struct BSubVertex : public BCell {
1532 typedef B_DCEL_MESH dcel_mesh_type;
1533 double icoor[2];
1534 };
1535
1536 struct SubVertex {
1537 GSubVertex g_subvertex;
1538 BSubVertex b_subvertex;
1539 SubVertex* g_previous_subvertex;
1540 SubVertex* g_next_subvertex;
1541 SubVertex* b_previous_subvertex;
1542 SubVertex* b_next_subvertex;
1543 };
1544
1545 struct CompSV {
1546 bool operator()(
1547 const std::pair<g_edge_type*, SubVertex*>& a,
1548 const std::pair<g_edge_type*, SubVertex*>& b) {
1549 if (a.first == b.first) {
1550 double ac;
1551 if (a.second->g_subvertex.dim == 1) {
1552 ac = a.second->g_subvertex.icoor[0];
1553 } else {
1554 assert(a.second->g_subvertex.dim == 0);
1555 if (a.second->g_subvertex.v == a.first->vertex) {
1556 ac = 1;
1557 } else {
1558 ac = 0;
1559 }
1560 }
1561 double bc;
1562 if (b.second->g_subvertex.dim == 1) {
1563 bc = b.second->g_subvertex.icoor[0];
1564 } else {
1565 assert(b.second->g_subvertex.dim == 0);
1566 if (b.second->g_subvertex.v == b.first->vertex) {
1567 bc = 1;
1568 } else {
1569 bc = 0;
1570 }
1571 }
1572 return ac < bc;
1573 } else {
1574 return a.first < b.first;
1575 }
1576 }
1577 };
1578
1579 void print_subvertex(const SubVertex& sv) const {
1580 std::cout << "sv=" << get_subvertex_id(&sv) << ", ptr = " << &sv << std::endl;
1581 const GSubVertex& gsv = sv.g_subvertex;
1582 std::cout << " g: dim = " << gsv.dim << ", id = ";
1583 switch (gsv.dim) {
1584 case 0:
1585 std::cout << gsv.v->id;
1586 break;
1587 case 1:
1588 std::cout << gsv.e->pair->vertex->id << ", " << gsv.e->vertex->id << ", coor = ("
1589 << gsv.icoor[0] << ")";
1590 break;
1591 case 2:
1592 if (gsv.f) {
1593 std::cout << gsv.f->id << ", coor = (" << gsv.icoor[0] << ", " << gsv.icoor[1]
1594 << ")";
1595 } else {
1596 std::cout << -1;
1597 }
1598 break;
1599 }
1600 std::cout << ", prev = " << get_subvertex_id(sv.g_previous_subvertex)
1601 << ", next = " << get_subvertex_id(sv.g_next_subvertex) << std::endl;
1602 const BSubVertex& bsv = sv.b_subvertex;
1603 std::cout << " b: dim = " << bsv.dim << ", id = ";
1604 switch (bsv.dim) {
1605 case 0:
1606 std::cout << bsv.v->id;
1607 break;
1608 case 1:
1609 std::cout << bsv.e->pair->vertex->id << ", " << bsv.e->vertex->id << ", coor = ("
1610 << bsv.icoor[0] << ")";
1611 break;
1612 case 2:
1613 if (bsv.f) {
1614 std::cout << bsv.f->id << ", coor = (" << bsv.icoor[0] << ", " << bsv.icoor[1]
1615 << ")";
1616 } else {
1617 std::cout << -1;
1618 }
1619 break;
1620 }
1621 std::cout << ", prev = " << get_subvertex_id(sv.b_previous_subvertex)
1622 << ", next = " << get_subvertex_id(sv.b_next_subvertex) << std::endl;
1623 }
1624
1625 void print_subvertex_gcoor(const SubVertex& sv) const {
1626 const GSubVertex& gsv = sv.g_subvertex;
1627 switch (gsv.dim) {
1628 case 0:
1629 std::cout << gsv.v->coor[0] << " " << gsv.v->coor[1] << " " << gsv.v->coor[2];
1630 break;
1631 case 1: {
1632 for (int i = 0; i != 3; ++i) {
1633 std::cout << gsv.e->vertex->coor[i] * gsv.icoor[0]
1634 + gsv.e->pair->vertex->coor[i] * (1 - gsv.icoor[0])
1635 << " ";
1636 }
1637 break;
1638 }
1639 case 2:
1640 if (gsv.f) {
1641 g_edge_type* e = gsv.f->edge;
1642 double* c0 = e->vertex->coor;
1643 e = e->next;
1644 double* c1 = e->vertex->coor;
1645 e = e->next;
1646 double* c2 = e->vertex->coor;
1647 const double h[3] = {
1648 1 - gsv.icoor[0] - gsv.icoor[1], gsv.icoor[0], gsv.icoor[1]};
1649 for (int i = 0; i != 3; ++i) {
1650 std::cout << c0[i] * h[0] + c1[i] * h[1] + c2[i] * h[2] << " ";
1651 }
1652 } else {
1653 std::cout << "-1 -1 -1";
1654 }
1655 break;
1656 }
1657 }
1658
1659 void print_subvertex_bcoor(const SubVertex& sv) const {
1660 const BSubVertex& gsv = sv.b_subvertex;
1661 switch (gsv.dim) {
1662 case 0:
1663 std::cout << gsv.v->coor[0] << " " << gsv.v->coor[1] << " " << gsv.v->coor[2];
1664 break;
1665 case 1: {
1666 for (int i = 0; i != 3; ++i) {
1667 std::cout << gsv.e->vertex->coor[i] * gsv.icoor[0]
1668 + gsv.e->pair->vertex->coor[i] * (1 - gsv.icoor[0])
1669 << " ";
1670 }
1671 break;
1672 }
1673 case 2:
1674 if (gsv.f) {
1675 b_edge_type* e = gsv.f->edge;
1676 double* c0 = e->vertex->coor;
1677 e = e->next;
1678 double* c1 = e->vertex->coor;
1679 e = e->next;
1680 double* c2 = e->vertex->coor;
1681 const double h[3] = {
1682 1 - gsv.icoor[0] - gsv.icoor[1], gsv.icoor[0], gsv.icoor[1]};
1683 for (int i = 0; i != 3; ++i) {
1684 std::cout << c0[i] * h[0] + c1[i] * h[1] + c2[i] * h[2] << " ";
1685 }
1686 } else {
1687 std::cout << "-1 -1 -1";
1688 }
1689 break;
1690 }
1691 }
1692
1693 void cluster_subvertex_on_dcel(SubVertex* a, SubVertex* b) {
1694 // blue mesh
1695 if (a->b_subvertex.dim > b->b_subvertex.dim) {
1696 if (a->b_subvertex.dim == 1) {
1697 if (a->b_subvertex.e->subvertex == a) {
1698 if (a->b_subvertex.e < a->b_subvertex.e->pair) {
1699 a->b_subvertex.e->subvertex =
1700 a->b_next_subvertex == b ? a : a->b_next_subvertex;
1701 } else {
1702 a->b_subvertex.e->subvertex =
1703 a->b_previous_subvertex == b ? a : a->b_previous_subvertex;
1704 }
1705 }
1706 if (a->b_subvertex.e->pair->subvertex == a) {
1707 if (a->b_subvertex.e->pair < a->b_subvertex.e) {
1708 a->b_subvertex.e->pair->subvertex =
1709 a->b_next_subvertex == b ? a : a->b_next_subvertex;
1710 } else {
1711 a->b_subvertex.e->pair->subvertex =
1712 a->b_previous_subvertex == b ? a : a->b_previous_subvertex;
1713 }
1714 }
1715 }
1716 if (a->b_previous_subvertex) {
1717 a->b_previous_subvertex->b_next_subvertex = a->b_next_subvertex;
1718 }
1719 if (a->b_next_subvertex) {
1720 a->b_next_subvertex->b_previous_subvertex = a->b_previous_subvertex;
1721 }
1722 if (b->b_subvertex.dim == 0) {
1723 if (b->b_subvertex.v->subvertex == b) { b->b_subvertex.v->subvertex = a; }
1724 b_edge_type* e = b->b_subvertex.v->edge;
1725 b_edge_type* e_end = e;
1726 do {
1727 if (e->subvertex) {
1728 SubVertex* sv = static_cast<SubVertex*>(e->subvertex);
1729 if (sv->b_next_subvertex == b) { sv->b_next_subvertex = a; }
1730 if (sv->b_previous_subvertex == b) { sv->b_previous_subvertex = a; }
1731 if (e->pair->subvertex == b) { e->pair->subvertex = a; }
1732 }
1733 e = e->pair->next;
1734 } while (e != e_end);
1735 } else if (b->b_subvertex.dim == 1) {
1736 if (b->b_subvertex.e->subvertex == b) { b->b_subvertex.e->subvertex = a; }
1737 if (b->b_subvertex.e->pair->subvertex == b) {
1738 b->b_subvertex.e->pair->subvertex = a;
1739 }
1740 }
1741 if (b->b_previous_subvertex && b->b_previous_subvertex->b_next_subvertex == b) {
1742 b->b_previous_subvertex->b_next_subvertex = a;
1743 }
1744 if (b->b_next_subvertex && b->b_next_subvertex->b_previous_subvertex == b) {
1745 b->b_next_subvertex->b_previous_subvertex = a;
1746 }
1747 std::swap(a->b_subvertex, b->b_subvertex);
1748 } else {
1749 if (b->b_subvertex.dim == 1) {
1750 if (b->b_subvertex.e->subvertex == b) {
1751 if (b->b_subvertex.e < b->b_subvertex.e->pair) {
1752 b->b_subvertex.e->subvertex = b->b_next_subvertex;
1753 } else {
1754 b->b_subvertex.e->subvertex = b->b_previous_subvertex;
1755 }
1756 }
1757 if (b->b_subvertex.e->pair->subvertex == b) {
1758 if (b->b_subvertex.e->pair < b->b_subvertex.e) {
1759 b->b_subvertex.e->pair->subvertex = b->b_next_subvertex;
1760 } else {
1761 b->b_subvertex.e->pair->subvertex = b->b_previous_subvertex;
1762 }
1763 }
1764 }
1765 if (b->b_previous_subvertex) {
1766 b->b_previous_subvertex->b_next_subvertex = b->b_next_subvertex;
1767 }
1768 if (b->b_next_subvertex) {
1769 b->b_next_subvertex->b_previous_subvertex = b->b_previous_subvertex;
1770 }
1771 if (b->b_subvertex.dim == 0) {
1772 if (b->b_subvertex.v->subvertex == b) { b->b_subvertex.v->subvertex = a; }
1773 b_edge_type* e = b->b_subvertex.v->edge;
1774 b_edge_type* e_end = e;
1775 do {
1776 if (e->subvertex) {
1777 SubVertex* sv = static_cast<SubVertex*>(e->subvertex);
1778 if (sv->b_next_subvertex == b) { sv->b_next_subvertex = a; }
1779 if (sv->b_previous_subvertex == b) { sv->b_previous_subvertex = a; }
1780 if (e->pair->subvertex == b) { e->pair->subvertex = a; }
1781 }
1782 e = e->pair->next;
1783 } while (e != e_end);
1784 }
1785 }
1786 b->b_next_subvertex = b->b_previous_subvertex = 0;
1787 b->b_subvertex.dim = -1;
1788
1789 // green mesh
1790 if (a->g_subvertex.dim > b->g_subvertex.dim) {
1791 if (a->g_subvertex.dim == 1) {
1792 if (a->g_subvertex.e->subvertex == a) {
1793 if (a->g_subvertex.e < a->g_subvertex.e->pair) {
1794 a->g_subvertex.e->subvertex =
1795 a->g_next_subvertex == b ? a : a->g_next_subvertex;
1796 } else {
1797 a->g_subvertex.e->subvertex =
1798 a->g_previous_subvertex == b ? a : a->g_previous_subvertex;
1799 }
1800 }
1801 if (a->g_subvertex.e->pair->subvertex == a) {
1802 if (a->g_subvertex.e->pair < a->g_subvertex.e) {
1803 a->g_subvertex.e->pair->subvertex =
1804 a->g_next_subvertex == b ? a : a->g_next_subvertex;
1805 } else {
1806 a->g_subvertex.e->pair->subvertex =
1807 a->g_previous_subvertex == b ? a : a->g_previous_subvertex;
1808 }
1809 }
1810 }
1811 if (a->g_previous_subvertex) {
1812 a->g_previous_subvertex->g_next_subvertex = a->g_next_subvertex;
1813 }
1814 if (a->g_next_subvertex) {
1815 a->g_next_subvertex->g_previous_subvertex = a->g_previous_subvertex;
1816 }
1817 if (b->g_subvertex.dim == 0) {
1818 if (b->g_subvertex.v->subvertex == b) { b->g_subvertex.v->subvertex = a; }
1819 g_edge_type* e = b->g_subvertex.v->edge;
1820 g_edge_type* e_end = e;
1821 do {
1822 if (e->subvertex) {
1823 SubVertex* sv = static_cast<SubVertex*>(e->subvertex);
1824 if (sv->g_next_subvertex == b) { sv->g_next_subvertex = a; }
1825 if (sv->g_previous_subvertex == b) { sv->g_previous_subvertex = a; }
1826 if (e->pair->subvertex == b) { e->pair->subvertex = a; }
1827 }
1828 e = e->pair->next;
1829 } while (e != e_end);
1830 } else if (b->g_subvertex.dim == 1) {
1831 if (b->g_subvertex.e->subvertex == b) { b->g_subvertex.e->subvertex = a; }
1832 if (b->g_subvertex.e->pair->subvertex == b) {
1833 b->g_subvertex.e->pair->subvertex = a;
1834 }
1835 }
1836 if (b->g_previous_subvertex && b->g_previous_subvertex->g_next_subvertex == b) {
1837 b->g_previous_subvertex->g_next_subvertex = a;
1838 }
1839 if (b->g_next_subvertex && b->g_next_subvertex->g_previous_subvertex == b) {
1840 b->g_next_subvertex->g_previous_subvertex = a;
1841 }
1842 std::swap(a->g_subvertex, b->g_subvertex);
1843 } else {
1844 if (b->g_subvertex.dim == 1) {
1845 if (b->g_subvertex.e->subvertex == b) {
1846 if (b->g_subvertex.e < b->g_subvertex.e->pair) {
1847 b->g_subvertex.e->subvertex = b->g_next_subvertex;
1848 } else {
1849 b->g_subvertex.e->subvertex = b->g_previous_subvertex;
1850 }
1851 }
1852 if (b->g_subvertex.e->pair->subvertex == b) {
1853 if (b->g_subvertex.e->pair < b->g_subvertex.e) {
1854 b->g_subvertex.e->pair->subvertex = b->g_next_subvertex;
1855 } else {
1856 b->g_subvertex.e->pair->subvertex = b->g_previous_subvertex;
1857 }
1858 }
1859 }
1860 if (b->g_previous_subvertex) {
1861 b->g_previous_subvertex->g_next_subvertex = b->g_next_subvertex;
1862 }
1863 if (b->g_next_subvertex) {
1864 b->g_next_subvertex->g_previous_subvertex = b->g_previous_subvertex;
1865 }
1866 if (b->g_subvertex.dim == 0) {
1867 if (b->g_subvertex.v->subvertex == b) { b->g_subvertex.v->subvertex = a; }
1868 g_edge_type* e = b->g_subvertex.v->edge;
1869 g_edge_type* e_end = e;
1870 do {
1871 if (e->subvertex) {
1872 SubVertex* sv = static_cast<SubVertex*>(e->subvertex);
1873 if (sv->g_next_subvertex == b) { sv->g_next_subvertex = a; }
1874 if (sv->g_previous_subvertex == b) { sv->g_previous_subvertex = a; }
1875 if (e->pair->subvertex == b) { e->pair->subvertex = a; }
1876 }
1877 e = e->pair->next;
1878 } while (e != e_end);
1879 }
1880 }
1881 b->g_next_subvertex = b->g_previous_subvertex = 0;
1882 b->g_subvertex.dim = -1;
1883 }
1884
1893 bool try_g_vertice_on_blue_face(
1894 b_face_type* bf, std::vector<void*> list_b_cell, g_edge_type* ge) {
1895 assert(bf != 0);
1896 std::stack<g_edge_type*> stack;
1897 stack.push(ge);
1898 bool res = false;
1899 while (!stack.empty()) {
1900 ge = stack.top();
1901 stack.pop();
1902 if (ge->vertex->subvertex) { continue; }
1903 BSubVertex bsv;
1904 double dist;
1905 if (g_point_to_b_face(*ge->vertex, *bf, bsv, dist)) {
1906 if (bsv.dim == 2) {
1907 SubVertex& rsv = new_sub_vertex();
1908 rsv.b_subvertex = bsv;
1909 rsv.g_subvertex.dim = 0;
1910 rsv.g_subvertex.v = ge->vertex;
1911 ge->vertex->subvertex = &rsv;
1912 } else {
1913 // Fix inconsistency introduced by the perturbation technique
1914 UnimplementedError() << THROW;
1915 }
1916 for (g_edge_type* e = ge->next; e->pair != ge; e = e->pair->next) {
1917 if (e->vertex->subvertex == 0) { stack.push(e); }
1918 }
1919 res = true;
1920 }
1921 }
1922 return res;
1923 }
1924
1925 bool intersect_triangle(
1926 const double orig[3], const double dir[3], const double vert0[3], const double vert1[3],
1927 const double vert2[3], double& u, double& v, double& t, double tol = -1) {
1928 if (tol < 0) { tol = tolerancing; }
1929 double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
1930 double det, inv_det;
1931 // find vectors for two edges sharing vert0
1932 sub_3(vert1, vert0, edge1);
1933 sub_3(vert2, vert0, edge2);
1934 // begin calculating determinant - also used to calculate U parameter
1935 outer_product_3(dir, edge2, pvec);
1936 // if determinant is near zero, ray lies in plane of triangle
1937 det = dot_3(edge1, pvec);
1938#ifdef TEST_CULL /* define TEST_CULL if culling is desired */
1939 if (det < epsilon) { return false; }
1940 /* calculate distance from vert0 to ray origin */
1941 sub_3(orig, vert0, tvec);
1942 /* calculate U parameter and test bounds */
1943 *u = dot_3(tvec, pvec);
1944 if (*u < -epsilon || *u > det + epsilon) { return false; }
1945 /* prepare to test V parameter */
1946 outer_product_3(tvec, edge1, qvec);
1947 *v = dot_3(dir, qvec);
1948 if (*v < -epsilon || *u + *v > det + epsilon) { return false; }
1949 /* calculate t, scale parameters, ray intersects triangle */
1950 *t = dot_3(edge2, qvec);
1951 inv_det = 1.0 / det;
1952 *t *= inv_det;
1953 *u *= inv_det;
1954 *v *= inv_det;
1955#else /* the non-culling branch */
1956 if (det > -tol && det < tol) { return false; }
1957 inv_det = 1.0 / det;
1958 /* calculate distance from vert0 to ray origin */
1959 sub_3(orig, vert0, tvec);
1960 /* calculate U parameter and test bounds */
1961 u = dot_3(tvec, pvec) * inv_det;
1962 if (u < -tol || u > 1.0 + tol) { return false; }
1963 /* prepare to test V parameter */
1964 outer_product_3(tvec, edge1, qvec);
1965 /* calculate V parameter and test bounds */
1966 v = dot_3(dir, qvec) * inv_det;
1967 if (v < -tol || u + v > 1.0 + tol) { return false; }
1968 /* calculate t, ray intersects triangle */
1969 t = dot_3(edge2, qvec) * inv_det;
1970#endif
1971 return true;
1972 }
1973
1974 bool g_point_to_b_face(
1975 g_vertex_type& g_vertex, b_face_type& face, BSubVertex& bsv, double& dist) {
1976 const double tolerancing = this->tolerancing * 0.01;
1977 const b_edge_type* vert0 = face.edge;
1978 const b_edge_type* vert1 = vert0->next;
1979 const b_edge_type* vert2 = vert1->next;
1980 if (!intersect_triangle(
1981 g_vertex.coor, g_vertex.normal, vert0->vertex->coor, vert1->vertex->coor,
1982 vert2->vertex->coor, bsv.icoor[0], bsv.icoor[1], dist, tolerancing)) {
1983 return false;
1984 }
1985
1986 const double x[3] = {bsv.icoor[0], bsv.icoor[1], 1 - x[0] - x[1]};
1987
1988 /*
1989 // a try to not collapse point with edge and vertex
1990 if (x[0] > -tolerancing && x[1] > -tolerancing && x[2] > -tolerancing) {
1991 if (bsv.icoor[0] < tolerancing)
1992 bsv.icoor[0] = tolerancing;
1993 if (bsv.icoor[0] > 1 - tolerancing)
1994 bsv.icoor[0] = 1 - tolerancing;
1995 if (bsv.icoor[1] < tolerancing)
1996 bsv.icoor[1] = tolerancing;
1997 if (bsv.icoor[1] > 1 - tolerancing)
1998 bsv.icoor[1] = 1 - tolerancing;
1999 if (bsv.icoor[0] + bsv.icoor[1] > 1 - tolerancing) {
2000 if (bsv.icoor[0] > bsv.icoor[1])
2001 bsv.icoor[0] = 1 - tolerancing - bsv.icoor[1];
2002 else
2003 bsv.icoor[1] = 1 - tolerancing - bsv.icoor[0];
2004 }
2005 bsv.dim = 2;
2006 bsv.f = &face;
2007 return true;
2008 }
2009 return false;
2010 */
2011
2012 if (x[0] > tolerancing && x[1] > tolerancing && x[2] > tolerancing) {
2013 bsv.dim = 2;
2014 bsv.f = &face;
2015 return true;
2016 }
2017 b_edge_type* e = face.edge;
2018 if (x[0] < tolerancing) {
2019 if (x[1] > tolerancing && x[2] > tolerancing) {
2020 bsv.dim = 1;
2021 bsv.e = e;
2022 bsv.icoor[0] = 1 - x[1];
2023 return true;
2024 } else if (x[1] < tolerancing) {
2025 bsv.dim = 0;
2026 bsv.v = e->vertex;
2027 return true;
2028 }
2029 }
2030 e = e->next;
2031 if (x[1] < tolerancing) {
2032 if (x[2] < tolerancing) {
2033 bsv.dim = 0;
2034 bsv.v = e->vertex;
2035 return true;
2036 } else {
2037 bsv.dim = 1;
2038 bsv.e = e;
2039 bsv.icoor[0] = x[0];
2040 return true;
2041 }
2042 }
2043 e = e->next;
2044 assert(x[2] < tolerancing);
2045 {
2046 if (x[0] < tolerancing) {
2047 bsv.dim = 0;
2048 bsv.v = e->vertex;
2049 } else {
2050 bsv.dim = 1;
2051 bsv.e = e;
2052 bsv.icoor[0] = x[1];
2053 }
2054 }
2055 return true;
2056 }
2057
2058 struct SystemSolveBVertexToGFaceProjection {
2059 typedef double value_type;
2060
2061 static const int size = 3;
2062
2063 void value(const double x[3], double y[3], double jacobien[9]) const {
2064 for (int i = 0; i != 3; ++i) {
2065 const double a1 = gp[0][i] - x[2] * gn[0][i];
2066 const double a2 = gp[1][i] - x[2] * gn[1][i];
2067 const double a3 = gp[2][i] - x[2] * gn[2][i];
2068 const double b = 1 - x[0] - x[1];
2069 y[i] = -bp[i] + b * a1 + x[0] * a2 + x[1] * a3;
2070 jacobien[i] = a2 - a1;
2071 jacobien[i + 3] = a3 - a1;
2072 jacobien[i + 6] = -b * gn[0][i] - x[0] * gn[1][i] - x[1] * gn[2][i];
2073 }
2074 }
2075 double bp[3];
2076 double gn[3][3];
2077 double gp[3][3];
2078 };
2079
2080 bool b_point_to_g_face(
2081 b_vertex_type& b_vertex, g_face_type& face, GSubVertex& gsv, double& dist) {
2082 SystemSolveBVertexToGFaceProjection sys;
2083 std::copy(b_vertex.coor, b_vertex.coor + 3, sys.bp);
2084 g_edge_type* e = face.edge;
2085 for (int i = 0; i != 3; ++i) {
2086 const g_vertex_type* v = e->vertex;
2087 std::copy(v->coor, v->coor + 3, sys.gp[i]);
2088 std::copy(v->normal, v->normal + 3, sys.gn[i]);
2089 e = e->next;
2090 }
2091 double x[3] = {};
2092 if (!solve_newton_raphson(sys, x)) { return false; }
2093 dist = x[2];
2094 x[2] = 1 - x[0] - x[1];
2095 if (x[0] < -tolerancing || x[1] < -tolerancing || x[2] < -tolerancing) { return false; }
2096 if (x[0] > tolerancing && x[1] > tolerancing && x[2] > tolerancing) {
2097 gsv.dim = 2;
2098 gsv.f = &face;
2099 gsv.icoor[0] = x[0];
2100 gsv.icoor[1] = x[1];
2101 return true;
2102 }
2103 e = face.edge;
2104 if (x[0] < tolerancing) {
2105 if (x[1] > tolerancing && x[2] > tolerancing) {
2106 gsv.dim = 1;
2107 gsv.e = e;
2108 gsv.icoor[0] = 1 - x[1];
2109 return true;
2110 } else if (x[1] < tolerancing) {
2111 gsv.dim = 0;
2112 gsv.v = e->vertex;
2113 return true;
2114 }
2115 }
2116 e = e->next;
2117 if (x[1] < tolerancing) {
2118 if (x[2] < tolerancing) {
2119 gsv.dim = 0;
2120 gsv.v = e->vertex;
2121 return true;
2122 } else {
2123 gsv.dim = 1;
2124 gsv.e = e;
2125 gsv.icoor[0] = x[0];
2126 return true;
2127 }
2128 }
2129 e = e->next;
2130 assert(x[2] < tolerancing);
2131 {
2132 if (x[0] < tolerancing) {
2133 gsv.dim = 0;
2134 gsv.v = e->vertex;
2135 } else {
2136 gsv.dim = 1;
2137 gsv.e = e;
2138 gsv.icoor[0] = x[1];
2139 }
2140 }
2141 return true;
2142 }
2143
2144 bool edge_intersection(
2145 b_edge_type& b_edge, g_edge_type& g_edge, SubVertex& sv, double& alpha, double& beta) {
2146 typedef long double T;
2147 const g_vertex_type& g_vertex0 = *g_edge.pair->vertex;
2148 const double* g0 = g_vertex0.coor;
2149 const double* d0 = g_vertex0.normal;
2150 const g_vertex_type& g_vertex1 = *g_edge.vertex;
2151 const double* g1 = g_vertex1.coor;
2152 const double* d1 = g_vertex1.normal;
2153 const double* b0 = b_edge.pair->vertex->coor;
2154 const double* b1 = b_edge.vertex->coor;
2155
2156 T b1_b0[3];
2157 sub_3(b1, b0, b1_b0);
2158 T d1_d0[3];
2159 sub_3(d1, d0, d1_d0);
2160 T g1_g0[3];
2161 sub_3(g1, g0, g1_g0);
2162 T b1_b0_cross_d1_d0[3];
2163 outer_product_3(b1_b0, d1_d0, b1_b0_cross_d1_d0);
2164 const T c2 = dot_3(b1_b0_cross_d1_d0, g1_g0);
2165 T b1_b0_cross_d0[3];
2166 outer_product_3(b1_b0, d0, b1_b0_cross_d0);
2167 T g0_b0[3];
2168 sub_3(g0, b0, g0_b0);
2169 const T c1 = dot_3(b1_b0_cross_d0, g1_g0) + dot_3(b1_b0_cross_d1_d0, g0_b0);
2170 const T c0 = dot_3(b1_b0_cross_d0, g0_b0);
2171
2172 int nb_beta;
2173 T l_beta[2];
2174 if (std::abs(c2) > epsilon) {
2175 T delta = c1 * c1 - 4 * c0 * c2;
2176 if (delta < -epsilon) {
2177 alpha = beta = std::numeric_limits<double>::quiet_NaN();
2178 // no solution -> no intersection
2179 return false;
2180 } else if (delta < epsilon) {
2181 delta = 0;
2182 }
2183 delta = std::sqrt(delta);
2184 T beta1 = (-c1 - delta) / (2 * c2);
2185 T beta2 = (-c1 + delta) / (2 * c2);
2186 if (beta1 < -tolerancing || beta1 > 1.0 + tolerancing) {
2187 nb_beta = 1;
2188 l_beta[0] = beta2;
2189 } else if (beta2 < -tolerancing || beta2 > 1.0 + tolerancing) {
2190 nb_beta = 1;
2191 l_beta[0] = beta1;
2192 } else {
2193 nb_beta = 2;
2194 l_beta[0] = beta1;
2195 l_beta[1] = beta2;
2196 }
2197 } else if (std::abs(c1) > epsilon) {
2198 nb_beta = 1;
2199 l_beta[0] = -c0 / c1;
2200 } else {
2201 // Infinite number of solution -> infinite number of intersection.
2202 // We set the intersection at the start or the end of the green edge
2203 nb_beta = 2;
2204 l_beta[0] = 0;
2205 l_beta[1] = 1;
2206 }
2207
2208 T l_alpha[2];
2209 bool l_alpha_nan[2];
2210 for (int j = 0; j != nb_beta; ++j) {
2211 T d0_beta_d1_d0[3];
2212 for (int i = 0; i != 3; ++i) { d0_beta_d1_d0[i] = d0[i] + l_beta[j] * d1_d0[i]; }
2213 T l[3];
2214 outer_product_3(g1_g0, d0_beta_d1_d0, l);
2215 const T l_b1_b0 = dot_3(l, b1_b0);
2216 if (std::abs(l_b1_b0) < epsilon) {
2217 l_alpha_nan[j] = true;
2218 } else {
2219 l_alpha_nan[j] = false;
2220 l_alpha[j] = dot_3(l, g0_b0) / l_b1_b0;
2221 }
2222 }
2223
2224 if (nb_beta == 2) {
2225 int i;
2226 if (l_alpha_nan[0] || l_alpha_nan[1]) {
2227 if (l_alpha_nan[0] && l_alpha_nan[1]) { return false; }
2228 i = l_alpha_nan[0] ? 1 : 0;
2229 } else if (l_alpha[0] > l_alpha[1]) {
2230 if (l_alpha[1] > -tolerancing) {
2231 i = 1;
2232 } else {
2233 i = 0;
2234 }
2235 } else if (l_alpha[0] > -tolerancing) {
2236 i = 0;
2237 } else {
2238 i = 1;
2239 }
2240 alpha = l_alpha[i];
2241 beta = l_beta[i];
2242 } else {
2243 beta = l_beta[0];
2244 if (l_alpha_nan[0]) { return false; }
2245 alpha = l_alpha[0];
2246 }
2247
2248 if (alpha < tolerancing) {
2249 if (alpha > -tolerancing) {
2250 sv.b_subvertex.dim = 0;
2251 if (b_edge.pair) {
2252 sv.b_subvertex.v = b_edge.pair->vertex;
2253 } else {
2254 b_edge_type* e = &b_edge;
2255 while (e->next != &b_edge) { e = e->next; }
2256 sv.b_subvertex.v = e->vertex;
2257 }
2258 } else {
2259 return false;
2260 }
2261 } else if (alpha > 1.0 - tolerancing) {
2262 if (alpha < 1.0 + tolerancing) {
2263 sv.b_subvertex.dim = 0;
2264 sv.b_subvertex.v = b_edge.vertex;
2265 } else {
2266 return false;
2267 }
2268 } else {
2269 sv.b_subvertex.dim = 1;
2270 sv.b_subvertex.e = &b_edge;
2271 sv.b_subvertex.icoor[0] = alpha;
2272 }
2273
2274 if (beta < tolerancing) {
2275 if (beta > -tolerancing) {
2276 sv.g_subvertex.dim = 0;
2277 sv.g_subvertex.v = g_edge.pair->vertex;
2278 } else {
2279 return false;
2280 }
2281 } else if (beta > 1.0 - tolerancing) {
2282 if (beta < 1.0 + tolerancing) {
2283 sv.g_subvertex.dim = 0;
2284 sv.g_subvertex.v = g_edge.vertex;
2285 } else {
2286 return false;
2287 }
2288 } else {
2289 sv.g_subvertex.dim = 1;
2290 sv.g_subvertex.e = &g_edge;
2291 sv.g_subvertex.icoor[0] = beta;
2292 }
2293
2294 return true;
2295 }
2296
2302 std::pair<b_vertex_type*, GSubVertex> green_parent_of_a_blue_vertex_by_brute_force() {
2303 typename B_DCEL_MESH::vertex_iterator bv = b_mesh->vertex_begin();
2304 g_face_type* g_face = 0;
2305 std::pair<b_vertex_type*, GSubVertex> res;
2306 while (g_face == 0 && bv != b_mesh->vertex_end()) {
2307 double g_face_dist = 1e100;
2308 for (typename G_DCEL_MESH::face_iterator gf = g_mesh->face_begin();
2309 gf != g_mesh->face_end(); ++gf) {
2310 GSubVertex gsv;
2311 double dist;
2312 if (b_point_to_g_face(*bv, *gf, gsv, dist)) {
2313 if (std::abs(dist) < g_face_dist) {
2314 g_face = &*gf;
2315 g_face_dist = std::abs(dist);
2316 res.second = gsv;
2317 res.first = &*bv;
2318 }
2319 }
2320 }
2321 ++bv;
2322 }
2323 if (g_face == 0) {
2324 res.first = &*b_mesh->vertex_begin();
2325 res.second.dim = 2;
2326 res.second.f = 0;
2327 std::fill_n(res.second.icoor, 2, 0);
2328 }
2329 return res;
2330 }
2331
2332 B_DCEL_MESH* b_mesh;
2333 G_DCEL_MESH* g_mesh;
2334 double tolerancing;
2335 double epsilon;
2336 b_edge_type* b_border_edge;
2337
2338 SubVertex& new_sub_vertex() {
2339 list_subvertex.push_back(SubVertex());
2340 SubVertex& res = list_subvertex.back();
2341 res.b_next_subvertex = 0;
2342 res.b_previous_subvertex = 0;
2343 res.g_next_subvertex = 0;
2344 res.g_previous_subvertex = 0;
2345 return res;
2346 }
2347
2348 size_t get_subvertex_id(const void* sv) const {
2349 for (size_t i = 0; i != list_subvertex.size(); ++i) {
2350 if (&list_subvertex[i] == sv) { return i; }
2351 }
2352 assert(sv == nullptr);
2353 return 11111111111111;
2354 }
2355
2356 typedef std::deque<SubVertex> list_subvertex_t;
2357 list_subvertex_t list_subvertex;
2358};
2359
2360template <typename DCEL_HALF_HEDGE, typename TRIANGLE>
2361void compute_overlaying_surface_meshes(
2362 const DCEL_HALF_HEDGE& b_mesh, const DCEL_HALF_HEDGE& g_mesh,
2363 std::vector<TRIANGLE>& overlaying_mesh) {
2364 ComputeOverlayingSurfaceMeshes<DCEL_HALF_HEDGE, TRIANGLE> tmp(b_mesh, g_mesh, overlaying_mesh);
2365 tmp.compute();
2366}
2367
2368} // namespace b2000
2369
2370#endif /* B2OVERLAYING_SURFACE_MESHES_H_ */
b2000::details::operator<
bool operator<(const csda< T1 > &a, const csda< T2 > &b)
Comparison of two csda numbers is performed on the real part only.
Definition b2csda.H:262
THROW
#define THROW
Definition b2exception.H:198
b2tensor_calculus.H
b2000::ComputeOverlayingSurfaceMeshes::OverlayingMeshIterator
Definition b2overlaying_surface_meshes.H:1054
b2000::ComputeOverlayingSurfaceMeshes
Definition b2overlaying_surface_meshes.H:121
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::sub_3
void sub_3(const T1 v1[3], const T2 v2[3], T3 dest[3])
Definition b2tensor_calculus.H:239
b2000::outer_product_3
void outer_product_3(const T1 a[3], const T2 b[3], T3 c[3])
Definition b2tensor_calculus.H:389
b2000::dot_3
T dot_3(const T a[3], const T b[3])
Definition b2tensor_calculus.H:328
b2000::UnimplementedError
GenericException< UnimplementedError_name > UnimplementedError
Definition b2exception.H:314