BRL-CAD
nmg_misc.c
Go to the documentation of this file.
1 /* N M G _ M I S C . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1993-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @addtogroup nmg */
21 /** @{ */
22 /** @file primitives/nmg/nmg_misc.c
23  *
24  * As the name implies, these are miscellaneous routines that work with
25  * the NMG structures.
26  *
27  */
28 /** @} */
29 
30 #include "common.h"
31 
32 #include <math.h>
33 #include <string.h>
34 #include "bio.h"
35 
36 #include "vmath.h"
37 #include "nmg.h"
38 #include "rtgeom.h"
39 #include "raytrace.h"
40 #include "nurb.h"
41 
42 #include "db.h" /* for debugging stuff at bottom */
43 
44 
45 int
46 nmg_snurb_calc_lu_uv_orient(const struct loopuse *lu)
47 {
48  struct edgeuse *eu;
49  int edge_count=0;
50  int edge_no;
51  vect_t area = VINIT_ZERO;
52  point_t *pts;
53 
54  NMG_CK_LOOPUSE(lu);
55 
56  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
57  bu_bomb("nmg_snurb_calc_lu_uv_orient: LU has no edges\n");
58 
59  if (*lu->up.magic_p != NMG_FACEUSE_MAGIC)
60  bu_bomb("nmg_snurb_calc_lu_uv_orient: LU is not part of a faceuse\n");
61 
62  NMG_CK_FACEUSE(lu->up.fu_p);
63  NMG_CK_FACE(lu->up.fu_p->f_p);
64 
65  if (*lu->up.fu_p->f_p->g.magic_p != NMG_FACE_G_SNURB_MAGIC)
66  bu_bomb("nmg_snurb_calc_lu_uv_orient: LU is not part of a SNURB face\n");
67 
68  /* count "pseudo-vertices" in loop */
69  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
70  struct edge_g_cnurb *eg;
71 
72  NMG_CK_EDGEUSE(eu);
73 
74  if (*eu->g.magic_p != NMG_EDGE_G_CNURB_MAGIC)
75  bu_bomb("nmg_snurb_calc_lu_uv_orient: EU on NURB face does not have edge_g_cnurb geometry\n");
76 
77  eg = eu->g.cnurb_p;
78  NMG_CK_EDGE_G_CNURB(eg);
79 
80  if (eg->order <= 0)
81  edge_count++;
82  else
83  edge_count += 5;
84  }
85 
86  /* allocate memory for "pseudo-vertices" */
87  pts = (point_t *)bu_calloc(edge_count, sizeof(point_t), "Orient_nurb_face_loops: pts");
88 
89  /* Assign uv geometry to each "pseudo-vertex" */
90  edge_no = 0;
91  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
92  struct edge_g_cnurb *eg;
93  struct vertexuse *vu;
94  struct vertexuse_a_cnurb *vg1;
95 
96  eg = eu->g.cnurb_p;
97 
98  if (eg->order <= 0) {
99  vu = eu->vu_p;
100  NMG_CK_VERTEXUSE(vu);
101  if (*vu->a.magic_p != NMG_VERTEXUSE_A_CNURB_MAGIC)
102  bu_bomb("Orient_nurb_face_loops: vertexuse in face_g_snurb faceuse doesn't have edge_g_cnurb attribute\n");
103  vg1 = vu->a.cnurb_p;
104  VMOVE(pts[edge_no], vg1->param);
105  edge_no++;
106  } else {
107  fastf_t t1, t2;
108  hpoint_t crv_pt;
109  int coords;
110  int i;
111 
112  t1 = eg->k.knots[0];
113  t2 = eg->k.knots[eg->k.k_size-1];
114  coords = RT_NURB_EXTRACT_COORDS(eg->pt_type);
115 
116  for (i = 0; coords > 0 && i < 5; i++) {
117  fastf_t t;
118 
119  t = t1 + (t2 - t1) * 0.2 * (fastf_t)i;
120 
121  VSETALLN(crv_pt, 0.0, coords);
122  rt_nurb_c_eval(eg, t, crv_pt);
123  if (RT_NURB_IS_PT_RATIONAL(eg->pt_type)) {
124  VSCALE(pts[edge_no], crv_pt, crv_pt[coords-1]);
125  } else {
126  VMOVE(pts[edge_no], crv_pt);
127  }
128  edge_no++;
129  }
130  }
131  }
132 
133  /* translate loop such that pts[0] is at (0, 0, 0) */
134  for (edge_no = 1; edge_no < edge_count; edge_no++) {
135  VSUB2(pts[edge_no], pts[edge_no], pts[0]);
136  pts[edge_no][Z] = 0.0;
137  }
138  VSETALL(pts[0], 0.0);
139 
140  /* calculate area of loop in uv-space */
141  for (edge_no = 1; edge_no < edge_count - 1; edge_no++) {
142  vect_t cross;
143 
144  VCROSS(cross, pts[edge_no], pts[edge_no+1]);
145  VADD2(area, area, cross);
146  }
147 
148  bu_free((char *)pts, "nmg_snurb_calc_lu_uv_orient: pts");
149 
150  if (area[Z] > 0.0)
151  return OT_SAME;
152  if (area[Z] < 0.0)
153  return OT_OPPOSITE;
154 
155  return OT_NONE;
156 }
157 
158 
159 void
160 nmg_snurb_fu_eval(const struct faceuse *fu, const fastf_t u, const fastf_t v, fastf_t *pt_on_srf)
161 {
162  struct face *f;
163  hpoint_t tmp_pt = HINIT_ZERO;
164 
165  NMG_CK_FACEUSE(fu);
166 
167  f = fu->f_p;
168  NMG_CK_FACE(f);
169  if (!f->g.magic_p) {
170  bu_log("nmg_snurb_fu_get_norm: face has no geometry (%p)\n", (void *)f);
171  bu_bomb("nmg_snurb_fu_get_norm: bad face\n");
172  }
173  if (*f->g.magic_p != NMG_FACE_G_SNURB_MAGIC) {
174  bu_log("nmg_snurb_fu_get_norm: face is not a NURB face (%p)\n", (void *)f);
175  bu_bomb("nmg_snurb_fu_get_norm: bad face\n");
176  }
177 
178  rt_nurb_s_eval(f->g.snurb_p, u, v, tmp_pt);
179 
180  if (RT_NURB_IS_PT_RATIONAL(f->g.snurb_p->pt_type)) {
181  double d;
182 
183  d = 1.0 / tmp_pt[W];
184  VSCALE(pt_on_srf, tmp_pt, d);
185  } else {
186  VMOVE(pt_on_srf, tmp_pt);
187  }
188 }
189 
190 
191 void
192 nmg_snurb_fu_get_norm(const struct faceuse *fu, const fastf_t u, const fastf_t v, fastf_t *norm)
193 {
194  struct face *f;
195 
196  NMG_CK_FACEUSE(fu);
197 
198  f = fu->f_p;
199  NMG_CK_FACE(f);
200  if (!f->g.magic_p) {
201  bu_log("nmg_snurb_fu_get_norm: face has no geometry (%p)\n", (void *)f);
202  bu_bomb("nmg_snurb_fu_get_norm: bad face\n");
203  }
204  if (*f->g.magic_p != NMG_FACE_G_SNURB_MAGIC) {
205  bu_log("nmg_snurb_fu_get_norm: face is not a NURB face (%p)\n", (void *)f);
206  bu_bomb("nmg_snurb_fu_get_norm: bad face\n");
207  }
208 
209  rt_nurb_s_norm(f->g.snurb_p, u, v, norm);
210 
211  if ((fu->orientation != OT_SAME) != (f->flip != 0))
212  VREVERSE(norm, norm);
213 }
214 
215 
216 void
217 nmg_snurb_fu_get_norm_at_vu(const struct faceuse *fu, const struct vertexuse *vu, fastf_t *norm)
218 {
219  struct vertexuse_a_cnurb *va;
220 
221  NMG_CK_FACEUSE(fu);
222  NMG_CK_VERTEXUSE(vu);
223 
224  if (!vu->a.magic_p) {
225  bu_log("nmg_snurb_fu_get_norm_at_vu: vertexuse does not have an attribute (%p)\n",
226  (void *)vu);
227  bu_bomb("nmg_snurb_fu_get_norm_at_vu: bad VU\n");
228  }
229 
230  if (*vu->a.magic_p != NMG_VERTEXUSE_A_CNURB_MAGIC) {
231  bu_log("nmg_snurb_fu_get_norm_at_vu: vertexuse does not have a cnurb attribute (%p)\n",
232  (void *)vu);
233  bu_bomb("nmg_snurb_fu_get_norm_at_vu: bad VU\n");
234  }
235 
236  va = vu->a.cnurb_p;
237  NMG_CK_VERTEXUSE_A_CNURB(va);
238 
239  nmg_snurb_fu_get_norm(fu, va->param[0], va->param[1], norm);
240 
241 }
242 
243 
244 void
245 nmg_find_zero_length_edges(const struct model *m)
246 {
247  struct bu_ptbl eu_tab;
248  struct edgeuse *eu;
249  int i;
250 
251  bu_ptbl_init(&eu_tab, 64, " &eu_tab");
252 
253  nmg_edgeuse_tabulate(&eu_tab, &m->magic);
254 
255  for (i=0; i<BU_PTBL_END(&eu_tab); i++) {
256  struct loopuse *lu;
257 
258  eu = (struct edgeuse *)BU_PTBL_GET(&eu_tab, i);
259  NMG_CK_EDGEUSE(eu);
260 
261  if (eu->vu_p->v_p != eu->eumate_p->vu_p->v_p)
262  continue;
263 
264  /* found a zero length edge */
265 
266  bu_log("Edgeuse %p (vp %p to vp %p)\n",
267  (void *)eu, (void *)eu->vu_p->v_p, (void *)eu->eumate_p->vu_p->v_p);
268  if (*eu->up.magic_p != NMG_LOOPUSE_MAGIC) {
269  bu_log("\tThis is a wire edge\n");
270  continue;
271  }
272 
273  lu = eu->up.lu_p;
274 
275  nmg_pr_lu_briefly(lu, "");
276  }
277 
278  bu_ptbl_free(&eu_tab);
279 }
280 
281 
282 /**
283  * Finds the topmost face in a shell (in given direction). Expects to
284  * have a translation table (variable "flags") for the model, and will
285  * ignore face structures that have their flag set in the table.
286  *
287  * dir must be X, Y, or Z
288  */
289 struct face *
290 nmg_find_top_face_in_dir(const struct shell *s, int dir, long int *flags)
291 {
292  fastf_t extreme_value=(-MAX_FASTF);
293  fastf_t extreme_slope=(-MAX_FASTF);
294  vect_t edge;
295  vect_t normal;
296  struct face *f_top=(struct face *)NULL;
297  struct edge *e_top=(struct edge *)NULL;
298  struct vertex *vp_top=(struct vertex *)NULL;
299  struct loopuse *lu;
300  struct faceuse *fu;
301  struct edgeuse *eu, *eu1;
302  struct vertexuse *vu;
303  struct vertex *v1, *v2;
304 
305  if (RTG.NMG_debug & DEBUG_BASIC)
306  bu_log("nmg_find_top_face_in_dir(s = %p, dir=%d, flags = %p)\n",
307  (void *)s, dir, (void *)flags);
308 
309  NMG_CK_SHELL(s);
310 
311  if (dir < X || dir > Z) {
312  bu_log("nmg_find_top_face_in_dir: illegal direction: %d\n", dir);
313  return (struct face *)NULL;
314  }
315 
316  /* find extreme vertex */
317  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
318  NMG_CK_FACEUSE(fu);
319 
320  /* skip flagged faceuses */
321  if (NMG_INDEX_TEST(flags, fu->f_p))
322  continue;
323  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
324  NMG_CK_LOOPUSE(lu);
325  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_EDGEUSE_MAGIC) {
326  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
327  NMG_CK_EDGEUSE(eu);
328 
329  if (eu->vu_p->v_p->vg_p->coord[dir] > extreme_value) {
330  extreme_value = eu->vu_p->v_p->vg_p->coord[dir];
331  vp_top = eu->vu_p->v_p;
332  }
333  }
334  }
335  }
336  }
337  if (vp_top == (struct vertex *)NULL) {
338  bu_log("Find_top_face_in_dir: Could not find extreme vertex");
339  return (struct face *)NULL;
340  }
341 
342  if (RTG.NMG_debug & DEBUG_BASIC)
343  bu_log("top vertex is %p at (%g %g %g)\n",
344  (void *)vp_top, V3ARGS(vp_top->vg_p->coord));
345 
346  /* find edge from vp_top with extreme slope in "dir" direction */
347  for (BU_LIST_FOR (vu, vertexuse, &vp_top->vu_hd)) {
348  struct vertexuse *vu1;
349 
350  NMG_CK_VERTEXUSE(vu);
351 
352  /* only consider edgeuses */
353  if (*vu->up.magic_p != NMG_EDGEUSE_MAGIC)
354  continue;
355 
356  eu = vu->up.eu_p;
357  NMG_CK_EDGEUSE(eu);
358 
359  if (RTG.NMG_debug & DEBUG_BASIC) {
360  bu_log("Checking edge (%g %g %g)<->(%g %g %g)\n",
361  V3ARGS(eu->vu_p->v_p->vg_p->coord),
362  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord));
363  }
364 
365  /* skip wire edges */
366  if (*eu->up.magic_p != NMG_LOOPUSE_MAGIC)
367  continue;
368 
369  /* skip wire loops */
370  if (*eu->up.lu_p->up.magic_p != NMG_FACEUSE_MAGIC)
371  continue;
372 
373  /* skip finished faces */
374  if (NMG_INDEX_TEST(flags, eu->up.lu_p->up.fu_p->f_p))
375  continue;
376 
377  /* skip edges from other shells */
378  if (nmg_find_s_of_eu(eu) != s)
379  continue;
380 
381  /* skip zero length edges */
382  if (eu->eumate_p->vu_p->v_p == vp_top)
383  continue;
384 
385  /* get vertex at other end of this edge */
386  vu1 = eu->eumate_p->vu_p;
387  NMG_CK_VERTEXUSE(vu1);
388 
389  /* make a unit vector in direction of edgeuse */
390  VSUB2(edge, vu1->v_p->vg_p->coord, vu->v_p->vg_p->coord);
391  VUNITIZE(edge);
392 
393  if (RTG.NMG_debug & DEBUG_BASIC) {
394  bu_log("Checking edge (%g %g %g)<->(%g %g %g)\n",
395  V3ARGS(eu->vu_p->v_p->vg_p->coord),
396  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord));
397  bu_log("\tedge direction = (%g %g %g)\n", V3ARGS(edge));
398  bu_log("\t\textreme slope = %g\n", extreme_slope);
399  }
400 
401  /* check against current maximum slope */
402  if (edge[dir] > extreme_slope) {
403  if (RTG.NMG_debug & DEBUG_BASIC)
404  bu_log("New top edge!\n");
405  extreme_slope = edge[dir];
406  e_top = eu->e_p;
407  }
408  }
409  if (e_top == (struct edge *)NULL) {
410  bu_log("Fix_normals: Could not find uppermost edge");
411  return (struct face *)NULL;
412  }
413 
414  eu = e_top->eu_p;
415 
416  v1 = eu->vu_p->v_p;
417  NMG_CK_VERTEX(v1);
418  v2 = eu->eumate_p->vu_p->v_p;
419  NMG_CK_VERTEX(v2);
420 
421  if (RTG.NMG_debug & DEBUG_BASIC)
422  bu_log("top EU is %p (%g %g %g) <-> (%g %g %g)\n",
423  (void *)eu, V3ARGS(v1->vg_p->coord),
424  V3ARGS(v2->vg_p->coord));
425 
426  /* now find the face containing edge between v1 nad v2
427  with "left-pointing vector" having the most extreme slope */
428  extreme_slope = (-MAX_FASTF);
429 
430  for (BU_LIST_FOR (vu, vertexuse, &v1->vu_hd)) {
431  vect_t left;
432  vect_t edge_dir;
433 
434  NMG_CK_VERTEXUSE(vu);
435  if (*vu->up.magic_p != NMG_EDGEUSE_MAGIC)
436  continue;
437 
438  eu1 = vu->up.eu_p;
439  NMG_CK_EDGEUSE(eu1);
440 
441  /* don't bother with anything but faces */
442  if (*eu1->up.magic_p != NMG_LOOPUSE_MAGIC)
443  continue;
444 
445  /* skip edges not between correct vertices */
446  if (eu1->eumate_p->vu_p->v_p != v2)
447  continue;
448 
449  lu = eu1->up.lu_p;
450  NMG_CK_LOOPUSE(lu);
451  if (*lu->up.magic_p != NMG_FACEUSE_MAGIC)
452  continue;
453 
454  /* fu is a faceuse containing "eu1" */
455  fu = lu->up.fu_p;
456  NMG_CK_FACEUSE(fu);
457 
458  /* skip faces from other shells and flagged faceuses */
459  if (fu->s_p != s || NMG_INDEX_TEST(flags, fu->f_p))
460  continue;
461 
462  /* make a vector in the direction of "eu1" */
463  if (RTG.NMG_debug & DEBUG_BASIC)
464  bu_log("test EU is %p (%g %g %g) <-> (%g %g %g)\n",
465  (void *)eu, V3ARGS(eu->vu_p->v_p->vg_p->coord),
466  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord));
467 
468  VSUB2(edge_dir, eu1->eumate_p->vu_p->v_p->vg_p->coord, eu1->vu_p->v_p->vg_p->coord);
469 
470  if (RTG.NMG_debug & DEBUG_BASIC)
471  bu_log("edge_dir is (%g %g %g)\n", V3ARGS(edge_dir));
472 
473  /* find the normal for this faceuse */
474  if (*fu->f_p->g.magic_p == NMG_FACE_G_PLANE_MAGIC) {
475  NMG_GET_FU_NORMAL(normal, fu);
476  } else if (*fu->f_p->g.magic_p == NMG_FACE_G_SNURB_MAGIC) {
477  nmg_snurb_fu_get_norm_at_vu(fu, eu1->vu_p, normal);
478  }
479 
480  if (RTG.NMG_debug & DEBUG_BASIC)
481  bu_log("fu normal is (%g %g %g)\n", V3ARGS(normal));
482 
483  /* normal cross edge direction gives vector in face */
484  VCROSS(left, normal, edge_dir);
485 
486  /* unitize to get slope */
487  VUNITIZE(left);
488  if (RTG.NMG_debug & DEBUG_BASIC) {
489  bu_log("left vector is (%g %g %g)\n", V3ARGS(left));
490  bu_log("\textreme slope in %d direction is %g\n", dir, extreme_slope);
491  }
492 
493  /* check against current most extreme slope */
494  if (left[dir] > extreme_slope) {
495  if (RTG.NMG_debug & DEBUG_BASIC)
496  bu_log("new f_top\n");
497  extreme_slope = left[dir];
498  f_top = fu->f_p;
499  }
500  }
501 
502  if (f_top == (struct face *)NULL) {
503  bu_log("Nmg_find_top_face_in_dir: Could not find uppermost face");
504  return (struct face *)NULL;
505  }
506 
507  if (RTG.NMG_debug & DEBUG_BASIC)
508  bu_log("nmg_find_top_face_in_dir: top face = %p, dir = %d, top vertex = %p (%g %g %g)\n",
509  (void *)f_top, dir,
510  (void *)vp_top, V3ARGS(vp_top->vg_p->coord));
511 
512  return f_top;
513 }
514 
515 
516 /**
517  * Finds the topmost face in a shell (in some direction). Expects to
518  * have a translation table (variable "flags") for the model, and will
519  * ignore face structures that have their flag set in the table.
520  *
521  * returns the top face in some direction.
522  *
523  * dir will be set to X, Y, or Z to indicate which top face was found.
524  */
525 struct face *
526 nmg_find_top_face(const struct shell *s, int *dir, long int *flags)
527 {
528  struct face *top_face;
529 
530  for (*dir=X; *dir<=Z; (*dir)++)
531  if ((top_face=nmg_find_top_face_in_dir(s, *dir, flags)) != (struct face *)NULL)
532  return top_face;
533 
534  /* give up!! */
535  bu_log("Nmg_find_top_face: Cannot find a top face\n");
536  *dir = (-32000); /* will hopefully cause an error if used */
537  return (struct face *)NULL;
538 
539 }
540 
541 
542 /**
543  * Passed an bu_ptbl structure containing one shell, this routine
544  * examines the other shells in the region to determine if any are
545  * void shells within the shell on the bu_ptbl list. Any such void
546  * shells found are added to the bu_ptbl list. The final result is a
547  * ptbl list of shells where the first shell on the list is the outer
548  * shell, and any additional shells one the list are void shells
549  * within that outer shell. This is a support routine for
550  * "nmg_find_outer_and_void_shells" and gets called for every outer
551  * shell in the region
552  */
553 struct top_face
554 {
555  struct shell *s;
556  struct face *f;
557  int dir;
558  vect_t normal;
559 };
560 
561 
562 HIDDEN void
563 nmg_assoc_void_shells(const struct nmgregion *r, struct bu_ptbl *shells, const struct bn_tol *ttol)
564 {
565  struct shell *outer_shell, *void_s, *s;
566  struct faceuse *fu;
567  struct loopuse *lu;
568  struct edgeuse *eu;
569  long *flags;
570  struct top_face *top_faces;
571  int total_shells=0;
572  int i;
573  int dir;
574 
575  NMG_CK_REGION(r);
576  BU_CK_PTBL(shells);
577  BN_CK_TOL(ttol);
578 
579  outer_shell = (struct shell *)BU_PTBL_GET(shells, 0);
580  NMG_CK_SHELL(outer_shell);
581 
582  /* count shells in region */
583  for (BU_LIST_FOR (s, shell, &r->s_hd))
584  total_shells++;
585 
586  /* make an array of shells and top faces */
587  top_faces = (struct top_face *)bu_calloc(total_shells, sizeof(struct top_face), "nmg_assoc_void_shells: top_faces");
588 
589  /* make flags array for use by "nmg_find_top_face" */
590  flags = (long *)bu_calloc(r->m_p->maxindex, sizeof(long), "nmg_find_outer_and_void_shells: flags");
591 
592  top_faces[0].s = outer_shell;
593  top_faces[0].f = nmg_find_top_face(outer_shell, &dir, flags);
594  top_faces[0].dir = dir;
595  fu = top_faces[0].f->fu_p;
596  if (fu->orientation != OT_SAME)
597  fu = fu->fumate_p;
598  NMG_GET_FU_NORMAL(top_faces[0].normal, fu);
599 
600  /* fill in top_faces array */
601  i = 0;
602  for (BU_LIST_FOR (s, shell, &r->s_hd)) {
603  if (s == outer_shell)
604  continue;
605 
606  top_faces[++i].s = s;
607  top_faces[i].f = nmg_find_top_face(s, &dir, flags);
608  top_faces[i].dir = dir;
609  if (top_faces[i].f == (struct face *)NULL)
610  bu_log("WARNING: nmg_assoc_void_shells() could not find top face for shell %p\n", (void *)s);
611  else {
612  fu = top_faces[i].f->fu_p;
613  if (fu->orientation != OT_SAME)
614  fu = fu->fumate_p;
615  NMG_GET_FU_NORMAL(top_faces[i].normal, fu);
616  }
617  }
618 
619  /* look for voids */
620  for (BU_LIST_FOR (void_s, shell, &r->s_hd)) {
621  struct face *void_f;
622  /* int wrong_void=0; */
623  vect_t normal;
624 
625  if (void_s == outer_shell)
626  continue;
627 
628  NMG_CK_SHELL(void_s);
629 
630  void_f = (struct face *)NULL;
631  for (i=0; i<total_shells; i++) {
632  if (top_faces[i].s == void_s) {
633  void_f = top_faces[i].f;
634  dir = top_faces[i].dir;
635  VMOVE(normal, top_faces[i].normal);
636  break;
637  }
638  }
639  if (void_f == (struct face *)NULL)
640  bu_bomb("nmg_assoc_void_shells: no top face for a shell\n");
641 
642  if (normal[dir] < 0.0) {
643  /* this is a void shell */
644  struct face *int_f;
645  struct shell *test_s;
646  int breakout=0;
647  int not_in_this_shell=0;
648 
649  /* this is a void shell
650  * but does it belong with outer_shell */
651  if (!V3RPP1_IN_RPP2(void_s->sa_p->min_pt, void_s->sa_p->max_pt, outer_shell->sa_p->min_pt, outer_shell->sa_p->max_pt)) {
652  continue;
653  }
654 
655  for (BU_LIST_FOR (fu, faceuse, &void_s->fu_hd)) {
656  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
657  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
658  continue;
659  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
660  int nmg_class;
661 
662  nmg_class = nmg_class_pt_s(eu->vu_p->v_p->vg_p->coord, outer_shell, 0, ttol);
663 
664  if (nmg_class == NMG_CLASS_AoutB) {
665  breakout = 1;
666  not_in_this_shell = 1;
667  break;
668  }
669  }
670  if (breakout)
671  break;
672  }
673  if (breakout)
674  break;
675  }
676 
677  if (not_in_this_shell)
678  continue;
679 
680  int_f = (struct face *)NULL;
681  for (i=0; i<total_shells; i++) {
682  if (top_faces[i].s == void_s) {
683  int_f = top_faces[i].f;
684  break;
685  }
686  }
687  if (int_f == (struct face *)NULL)
688  bu_bomb("nmg_assoc_void_shells: no top face for a shell\n");
689 
690  /* Make sure there are no other external shells between these two */
691  for (BU_LIST_FOR (test_s, shell, &r->s_hd)) {
692  vect_t test_norm;
693  struct face *test_f;
694  int test_dir = 0;
695 
696  /* don't check against the outer shell or the candidate void shell */
697  if (test_s == void_s || test_s == outer_shell)
698  continue;
699 
700  /* find top face for the test shell */
701  test_f = (struct face *)NULL;
702  for (i=0; i<total_shells; i++) {
703  if (top_faces[i].s == test_s) {
704  test_f = top_faces[i].f;
705  test_dir = top_faces[i].dir;
706  VMOVE(test_norm, top_faces[i].normal);
707  break;
708  }
709  }
710  if (test_f == (struct face *)NULL)
711  bu_bomb("nmg_assoc_void_shells: no top face for a shell\n");
712 
713  /* skip test shells that are void shells */
714  if (test_norm[test_dir] < 0.0)
715  continue;
716 
717  /* if the void shell is not within the test shell, continue */
718  if (!V3RPP1_IN_RPP2(void_s->sa_p->min_pt, void_s->sa_p->max_pt, test_s->sa_p->min_pt, test_s->sa_p->max_pt))
719  continue;
720 
721  /* the void shell may be within this shell */
722  /* XXXX Need code here to check if candidate void shell (void_s)
723  * XXXX is within test shell (test_s) and test shell is
724  * XXXX is within outer shell (outer_shell)
725  if (void_s in test_s and test_s in outer_shell) {
726  wrong_void = 1;
727  break;
728  }
729  */
730  }
731  /*
732  if (wrong_void) {
733  continue;
734  }
735  */
736 
737  /* This void shell belongs with shell outer_s
738  * add it to the list of shells */
739  bu_ptbl_ins(shells, (long *)void_s);
740  }
741  }
742  bu_free((char *)flags, "nmg_assoc_void_shells: flags");
743 }
744 
745 
746 /**
747  * This routine takes a region and constructs an array of bu_ptbl
748  * lists. A list is created for each outer shell, and that shell is
749  * the first item on the list. Additional shells on any list are void
750  * shells within that lists outer shell. This routine calls
751  * "nmg_decompose_shell" for every shell in the region, so the
752  * original region topology may be changed to accomplish this. No
753  * geometry is altered.
754  */
755 int
756 nmg_find_outer_and_void_shells(struct nmgregion *r, struct bu_ptbl ***shells, const struct bn_tol *tol)
757 {
758  struct bu_ptbl *outer_shells;
759  struct shell *s;
760  int i;
761  int total_shells=0;
762  int outer_shell_count;
763  int re_bound=0;
764  int dir;
765  long *flags;
766 
767  NMG_CK_REGION(r);
768  BN_CK_TOL(tol);
769 
770  /* Decompose shells */
771  outer_shells = (struct bu_ptbl *)bu_malloc(sizeof(struct bu_ptbl), "nmg_find_outer_and_void_shells: outer_shells");
772  bu_ptbl_init(outer_shells, 64, " outer_shells ");
773  for (BU_LIST_FOR (s, shell, &r->s_hd)) {
774  NMG_CK_SHELL(s);
775  bu_ptbl_ins(outer_shells, (long *)s);
776  }
777  for (i=0; i<BU_PTBL_END(outer_shells); i++) {
778  s = (struct shell *)BU_PTBL_GET(outer_shells, i);
779  if (nmg_decompose_shell(s, tol) > 1)
780  re_bound = 1;
781  }
782  bu_ptbl_reset(outer_shells);
783 
784  if (re_bound)
785  nmg_region_a(r, tol);
786 
787  for (BU_LIST_FOR (s, shell, &r->s_hd))
788  total_shells++;
789 
790  flags = (long *)bu_calloc(r->m_p->maxindex, sizeof(long), "nmg_find_outer_and_void_shells: flags");
791 
792  for (BU_LIST_FOR (s, shell, &r->s_hd)) {
793  struct face *f;
794  struct faceuse *fu;
795  vect_t normal;
796 
797  f = (struct face *)NULL;
798  for (dir = X; dir <= Z; dir++) {
799  if ((f = nmg_find_top_face_in_dir(s, dir, flags)) == (struct face *)NULL)
800  continue;
801 
802  fu = f->fu_p;
803  if (fu->orientation != OT_SAME)
804  fu = fu->fumate_p;
805  if (fu->orientation != OT_SAME)
806  bu_bomb("nmg_find_outer_and_void_shells: Neither faceuse nor mate have OT_SAME orient\n");
807 
808  NMG_GET_FU_NORMAL(normal, fu);
809  if (normal[dir] >= 0.0) {
810  bu_ptbl_ins(outer_shells, (long *)s); /* outer shell */
811  break;
812  }
813  }
814 
815  if (f == (struct face *)NULL) {
816  bu_bomb("nmg_find_outer_and_void_shells: cannot find top face in a shell\n");
817  }
818  }
819 
820  /* outer_shells is now a list of all the outer shells in the region */
821  outer_shell_count = BU_PTBL_END(outer_shells);
822 
823  *shells = (struct bu_ptbl **)bu_calloc(BU_PTBL_END(outer_shells), sizeof(struct bu_ptbl *) ,
824  "nmg_find_outer_and_void_shells: shells");
825  for (i=0; i<BU_PTBL_END(outer_shells); i++) {
826  BU_ALLOC((*shells)[i], struct bu_ptbl);
827 
828  bu_ptbl_init((*shells)[i], 64, "(*shells)[i]");
829  BU_CK_PTBL((*shells)[i]);
830  bu_ptbl_ins((*shells)[i], BU_PTBL_GET(outer_shells, i));
831  if (outer_shell_count != total_shells) /* must be some void shells */
832  nmg_assoc_void_shells(r, (*shells)[i], tol);
833  }
834 
835  bu_free((char *)flags, "nmg_find_outer_and_void_shells: flags");
836  bu_ptbl_free(outer_shells);
837  return outer_shell_count;
838 }
839 
840 
841 /**
842  * Sets the "is_real" flag on all edges at or below the pointer
843  * passed. Returns the number of flags set.
844  */
845 int
846 nmg_mark_edges_real(const uint32_t *magic_p)
847 {
848  struct bu_ptbl edges;
849  int i, count;
850 
851  nmg_edge_tabulate(&edges, magic_p);
852 
853  count = BU_PTBL_END(&edges);
854  for (i=0; i<count; i++) {
855  struct edge *e;
856 
857  e = (struct edge *)BU_PTBL_GET(&edges, i);
858  NMG_CK_EDGE(e);
859 
860  e->is_real = 1;
861  }
862 
863  bu_ptbl_free(&edges);
864 
865  return count;
866 }
867 
868 
869 /**
870  * Tabulates all vertices in faces that use fg
871  */
872 void
873 nmg_tabulate_face_g_verts(struct bu_ptbl *tab, const struct face_g_plane *fg)
874 {
875  struct face *f;
876 
877  NMG_CK_FACE_G_PLANE(fg);
878 
879  bu_ptbl_init(tab, 64, " tab");
880 
881  /* loop through all faces using fg */
882  for (BU_LIST_FOR (f, face, &fg->f_hd)) {
883  struct faceuse *fu;
884  struct loopuse *lu;
885 
886  NMG_CK_FACE(f);
887 
888  /* get one of the two uses of this face */
889  fu = f->fu_p;
890  NMG_CK_FACEUSE(fu);
891 
892  /* Visit each loop in this faceuse */
893  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
894  NMG_CK_LOOPUSE(lu);
895 
896  /* include loops of a single vertex */
897  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
898  struct vertexuse *vu;
899  struct vertex *v;
900 
901  vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
902  NMG_CK_VERTEXUSE(vu);
903  v = vu->v_p;
904  NMG_CK_VERTEX(v);
905 
906  /* insert vertex into table */
907  bu_ptbl_ins_unique(tab, (long *)v);
908  } else {
909  struct edgeuse *eu;
910 
911  /* visit each edgeuse in the loop */
912  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
913  struct vertexuse *vu;
914  struct vertex *v;
915 
916  NMG_CK_EDGEUSE(eu);
917  vu = eu->vu_p;
918  NMG_CK_VERTEXUSE(vu);
919  v = vu->v_p;
920  NMG_CK_VERTEX(v);
921 
922  /* insert vertex into table */
923  bu_ptbl_ins_unique(tab, (long *)v);
924  }
925  }
926  }
927  }
928 }
929 
930 
931 /**
932  * Intersects all faces in a shell with all other faces in the same
933  * shell Intended for use after extrusion
934  */
935 void
936 nmg_isect_shell_self(struct shell *s, const struct bn_tol *tol)
937 {
938  struct model *m;
939  struct nmgregion *r;
940  struct shell *s_fu;
941  struct faceuse *fu;
942  struct bu_ptbl fus;
943  int fu_no;
944  int fu2_no;
945 
946  NMG_CK_SHELL(s);
947  BN_CK_TOL(tol);
948 
949  m = nmg_find_model(&s->l.magic);
950  NMG_CK_MODEL(m);
951 
952  nmg_vmodel(m);
953 
954  r = s->r_p;
955  NMG_CK_REGION(r);
956 
957  s_fu = nmg_msv(r);
958  NMG_CK_SHELL(s_fu);
959 
960  bu_ptbl_init(&fus, 64, " &fus ");
961 
962  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
963  NMG_CK_FACEUSE(fu);
964 
965  if (fu->orientation == OT_SAME)
966  bu_ptbl_ins(&fus, (long *)fu);
967  }
968 
969  /* intersect each face with every other face in the shell */
970  for (fu_no=0; fu_no < BU_PTBL_END(&fus); fu_no ++) {
971  struct faceuse *fu2;
972 
973  fu = (struct faceuse *)BU_PTBL_GET(&fus, fu_no);
974 
975  NMG_CK_FACEUSE(fu);
976 
977  /* move fu to another shell to avoid radial edge problems */
978  nmg_mv_fu_between_shells(s_fu, s, fu);
979 
980  /* consider intersection this faceuse with all the faceuses
981  * after it in the list
982  */
983  for (fu2_no=fu_no+1; fu2_no < BU_PTBL_END(&fus); fu2_no++) {
984  struct face *f, *f2;
985 
986  fu2 = (struct faceuse *)BU_PTBL_GET(&fus, fu2_no);
987 
988  if (RTG.NMG_debug & DEBUG_BASIC)
989  bu_log("nmg_extrude_cleanup: fu=%p, fu2=%p\n", (void *)fu, (void *)fu2);
990 
991  /* skip faceuses radial to fu or not OT_SAME */
992  if (fu2->orientation != OT_SAME || nmg_faces_are_radial(fu, fu2))
993  continue;
994 
995  f = fu->f_p;
996  f2 = fu2->f_p;
997 
998  /* skip faceuse pairs that don't have overlapping BB's */
999  if (!V3RPP_OVERLAP(f->min_pt, f->max_pt, f2->min_pt, f2->max_pt))
1000  continue;
1001 
1002  if (RTG.NMG_debug & DEBUG_BASIC)
1003  bu_log("nmg_extrude_cleanup: calling nmg_isect_two_generic_faces(fu=%p, fu2=%p)\n",
1004  (void *)fu, (void *)fu2);
1005 
1006  nmg_isect_two_generic_faces(fu, fu2, tol);
1007  }
1008  /* move fu back where it belongs */
1009  while (BU_LIST_NON_EMPTY(&s_fu->fu_hd)) {
1010  struct faceuse *fu_tmp;
1011 
1012  fu_tmp = BU_LIST_FIRST(faceuse, &s_fu->fu_hd);
1013  NMG_CK_FACEUSE(fu_tmp);
1014 
1015  if (RTG.NMG_debug & DEBUG_BASIC)
1016  bu_log("nmg_extrude_cleanup: moving fu %p back\n", (void *)fu_tmp);
1017 
1018  nmg_mv_fu_between_shells(s, s_fu, fu_tmp);
1019  }
1020  }
1021 
1022  /* get rid of the temporary shell */
1023  nmg_ks(s_fu);
1024 
1025  bu_ptbl_free(&fus);
1026 }
1027 
1028 
1029 /**
1030  * Traverse radial edgeuse around specified edgeuse looking for one
1031  * that meets optional restrictions. If a shell is specified only
1032  * edgeuse from that shell will be considered. If wires is non-zero,
1033  * wire edges will be considered, otherwise, wire edges are ignored.
1034  *
1035  * returns:
1036  * radial edgeuse
1037  */
1038 struct edgeuse *
1039 nmg_next_radial_eu(const struct edgeuse *eu, const struct shell *s, const int wires)
1040 {
1041  struct edgeuse *ret_eu;
1042 
1043  NMG_CK_EDGEUSE(eu);
1044  if (s)
1045  NMG_CK_SHELL(s);
1046 
1047  if (s && nmg_find_s_of_eu(eu) != s)
1048  bu_bomb("nmg_find_radial_eu: eu is not in specified shell\n");
1049 
1050  if (!wires && !nmg_find_fu_of_eu(eu))
1051  bu_bomb("nmg_find_radial_eu: wire edges not specified, but eu is a wire!!\n");
1052 
1053  ret_eu = eu->eumate_p->radial_p;
1054  while (
1055  (!wires & (nmg_find_fu_of_eu(ret_eu) == (struct faceuse *)NULL))
1056  ||
1057  ((s != (struct shell *)NULL) &&
1058  nmg_find_s_of_eu(ret_eu) != s)
1059  )
1060  ret_eu = ret_eu->eumate_p->radial_p;
1061 
1062  return ret_eu;
1063 }
1064 
1065 
1066 /**
1067  * Traverse radial edgeuse around specified edgeuse in opposite
1068  * direction from nmg_next_radial_eu, looking for one that meets
1069  * optional restrictions. If a shell is specified only edgeuse from
1070  * that shell will be considered. If wires is non-zero, wire edges
1071  * will be considered, otherwise, wire edges are ignored.
1072  *
1073  * returns:
1074  * radial edgeuse
1075  */
1076 struct edgeuse *
1077 nmg_prev_radial_eu(const struct edgeuse *eu, const struct shell *s, const int wires)
1078 {
1079  struct edgeuse *ret_eu;
1080 
1081  NMG_CK_EDGEUSE(eu);
1082  if (s)
1083  NMG_CK_SHELL(s);
1084 
1085  if (s && nmg_find_s_of_eu(eu) != s)
1086  bu_bomb("nmg_find_radial_eu: eu is not in specified shell\n");
1087 
1088  if (!wires && !nmg_find_fu_of_eu(eu))
1089  bu_bomb("nmg_find_radial_eu: wire edges not specified, but eu is a wire!!\n");
1090 
1091  ret_eu = eu->radial_p->eumate_p;
1092  while ((!wires & (nmg_find_fu_of_eu(ret_eu) == (struct faceuse *)NULL)) ||
1093  ((s != (struct shell *)NULL) && nmg_find_s_of_eu(ret_eu) != s))
1094  ret_eu = ret_eu->radial_p->eumate_p;
1095 
1096  return ret_eu;
1097 }
1098 
1099 
1100 /**
1101  * Counts the number of faces (actually, the number of radial edgeuse/mate pairs)
1102  * around eu. If s is specified, only edgeuses in shell s are counted. Wire
1103  * edgeuses are not counted.
1104  *
1105  * returns:
1106  * number of edgeuse/mate pairs radially around eu that meet restrictions
1107  */
1108 int
1109 nmg_radial_face_count(const struct edgeuse *eu, const struct shell *s)
1110 {
1111  int face_count=1;
1112  struct edgeuse *eu1;
1113 
1114  NMG_CK_EDGEUSE(eu);
1115  if (s)
1116  NMG_CK_SHELL(s);
1117 
1118  /* count radial faces on this edge */
1119  eu1 = eu->eumate_p->radial_p;
1120  while (eu1 != eu && eu1 != eu->eumate_p) {
1121  /* ignore other shells and don't count wires */
1122  if ((!s || nmg_find_s_of_eu(eu1) == s) &&
1123  nmg_find_fu_of_eu(eu1) != (struct faceuse *)NULL)
1124  face_count++;
1125  eu1 = eu1->eumate_p->radial_p;
1126  }
1127 
1128  return face_count;
1129 }
1130 
1131 
1132 /**
1133  * Looks at every eu in OT_SAME fu's. If any eu has no radials, then
1134  * it must be the edge of a dangling face and therefore the edge of an
1135  * opening.
1136  *
1137  * returns:
1138  * 0 - O.K.
1139  * 1 - found a hole
1140  */
1141 int
1142 nmg_check_closed_shell(const struct shell *s, const struct bn_tol *tol)
1143 {
1144  struct faceuse *fu;
1145 
1146  NMG_CK_SHELL(s);
1147  BN_CK_TOL(tol);
1148 
1149  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
1150  struct loopuse *lu;
1151 
1152  NMG_CK_FACEUSE(fu);
1153 
1154  if (fu->orientation != OT_SAME)
1155  continue;
1156 
1157  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
1158  struct edgeuse *eu;
1159 
1160  NMG_CK_LOOPUSE(lu);
1161 
1162  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
1163  continue;
1164 
1165  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
1166  struct edgeuse *next_eu;
1167 
1168  next_eu = nmg_next_radial_eu(eu, s, 0);
1169  if (next_eu == eu || next_eu == eu->eumate_p)
1170  return 1;
1171  }
1172  }
1173  }
1174 
1175  return 0;
1176 }
1177 
1178 
1179 /**
1180  * Moves lu from src faceuse to dest faceuse
1181  *
1182  * returns:
1183  * 0 - All is well
1184  * 1 - src faceuse is now empty
1185  */
1186 int
1187 nmg_move_lu_between_fus(struct faceuse *dest, struct faceuse *src, struct loopuse *lu)
1188 {
1189  struct loopuse *lumate;
1190  int src_is_empty;
1191 
1192  NMG_CK_FACEUSE(dest);
1193  NMG_CK_FACEUSE(dest->fumate_p);
1194  NMG_CK_FACEUSE(src);
1195  NMG_CK_FACEUSE(src->fumate_p);
1196  NMG_CK_LOOPUSE(lu);
1197 
1198  if (RTG.NMG_debug & DEBUG_BASIC)
1199  bu_log("nmg_move_lu_between_fus(dest=%p, src=%p, lu=%p)\n",
1200  (void *)dest, (void *)src, (void *)lu);
1201 
1202  if (lu->up.fu_p != src) {
1203  bu_log("nmg_move_lu_between_fus(dest=%p, src=%p, lu=%p)\n",
1204  (void *)dest, (void *)src, (void *)lu);
1205  bu_bomb("\tlu is not in src faceuse\n");
1206  }
1207 
1208  if (dest == src)
1209  return 0;
1210 
1211  lumate = lu->lumate_p;
1212  NMG_CK_LOOPUSE(lumate);
1213 
1214  /* remove lu from src faceuse */
1215  BU_LIST_DEQUEUE(&lu->l);
1216  src_is_empty = BU_LIST_IS_EMPTY(&src->lu_hd);
1217 
1218  /* remove lumate from src faceuse mate */
1219  BU_LIST_DEQUEUE(&lumate->l);
1220  if (src_is_empty != BU_LIST_IS_EMPTY(&src->fumate_p->lu_hd)) {
1221  bu_log("nmg_move_lu_between_fus(dest=%p, src=%p, lu=%p)\n",
1222  (void *)dest, (void *)src, (void *)lu);
1223  if (src_is_empty)
1224  bu_bomb("\tsrc faceuse contains only lu, but src->fumate_p has more!!\n");
1225  bu_bomb("\tsrc->fumate_p faceuse contains only lu->lumate_p, but src has more!!\n");
1226  }
1227 
1228  /* add lu to dest faceuse */
1229  BU_LIST_INSERT(&dest->lu_hd, &lu->l);
1230 
1231  /* add lumate to dest mate */
1232  BU_LIST_INSERT(&dest->fumate_p->lu_hd, &lumate->l);
1233 
1234  /* adjust up pointers */
1235  lu->up.fu_p = dest;
1236  lumate->up.fu_p = dest->fumate_p;
1237 
1238  return src_is_empty;
1239 }
1240 
1241 
1242 /**
1243  * Calculate the plane equation of a loop using Newell's Method (See
1244  * "Graphics Gems III", David Kirk editor, Academic Press, Inc. 1992)
1245  *
1246  * If the loop orientation is OT_OPPOSITE, the normal of the plane is
1247  * reversed.
1248  */
1249 void
1250 nmg_loop_plane_newell(const struct loopuse *lu, fastf_t *pl)
1251 {
1252  struct edgeuse *eu;
1253  double hmin, hmax;
1254  plane_t pl_tmp = HINIT_ZERO;
1255 
1256  NMG_CK_LOOPUSE(lu);
1257 
1258  HSETALL(pl, 0.0);
1259 
1260  /* make sure we have a loop of edges */
1261  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
1262  return;
1263 
1264  /* check if this loop is a crack */
1265  if (nmg_loop_is_a_crack(lu))
1266  return;
1267 
1268  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
1269  return;
1270 
1271  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
1272  struct edgeuse *eu_next;
1273  struct vertex_g *vg;
1274  struct vertex_g *vg_next;
1275 
1276  vg = eu->vu_p->v_p->vg_p;
1277  eu_next = BU_LIST_PNEXT_CIRC(edgeuse, &eu->l);
1278  vg_next = eu_next->vu_p->v_p->vg_p;
1279 
1280  pl_tmp[X] += (vg->coord[Y] - vg_next->coord[Y]) * (vg->coord[Z] + vg_next->coord[Z]);
1281  pl_tmp[Y] += (vg->coord[Z] - vg_next->coord[Z]) * (vg->coord[X] + vg_next->coord[X]);
1282  pl_tmp[Z] += (vg->coord[X] - vg_next->coord[X]) * (vg->coord[Y] + vg_next->coord[Y]);
1283  }
1284 
1285  VUNITIZE(pl_tmp);
1286  VMOVE(pl, pl_tmp);
1287 
1288  hmin = MAX_FASTF;
1289  hmax = (-hmin);
1290 
1291  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
1292  struct vertex_g *vg;
1293  fastf_t htest;
1294 
1295  vg = eu->vu_p->v_p->vg_p;
1296  htest = VDOT(vg->coord, pl);
1297  if (htest > hmax)
1298  hmax = htest;
1299  if (htest < hmin)
1300  hmin = htest;
1301  }
1302 
1303  pl[H] = (hmax + hmin)/2.0;
1304 
1305  if (lu->orientation == OT_OPPOSITE)
1306  HREVERSE(pl, pl);
1307 }
1308 
1309 
1310 /**
1311  * Calculates a plane equation and the area of a loop
1312  *
1313  * returns:
1314  * the area of the loop
1315  * less than zero indicates an error
1316  *
1317  * pl is assigned the plane equation for the loop
1318  */
1319 fastf_t
1320 nmg_loop_plane_area(const struct loopuse *lu, fastf_t *pl)
1321 {
1322  fastf_t area;
1323  fastf_t pt_count=0.0;
1324  fastf_t pt_dot_plane=0.0;
1325  plane_t plane = HINIT_ZERO;
1326  struct edgeuse *eu;
1327  vect_t trans;
1328 
1329  NMG_CK_LOOPUSE(lu);
1330 
1331  /* make sure we have a loop of edges */
1332  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
1333  return (fastf_t)(-1.0);
1334 
1335  /* check if this loop is a crack */
1336  if (nmg_loop_is_a_crack(lu))
1337  return (fastf_t)(-1.0);
1338 
1339  /* calculate a translation to put one vertex at the origin
1340  * not necessary, but good for accuracy.
1341  * Also, origin must be in plane of loop for this
1342  * method to work
1343  */
1344  eu = BU_LIST_FIRST(edgeuse, &lu->down_hd);
1345  NMG_CK_VERTEXUSE(eu->vu_p);
1346  NMG_CK_VERTEX(eu->vu_p->v_p);
1347  NMG_CK_VERTEX_G(eu->vu_p->v_p->vg_p);
1348  VMOVE(trans, eu->vu_p->v_p->vg_p->coord);
1349 
1350  /* Calculate area and plane normal.
1351  * Cross product of each pair of vertices gives twice
1352  * the area of the triangle formed by the origin and
1353  * the two vertices. (positive if counter-clockwise,
1354  * negative if clockwise). In counter_clockwise case,
1355  * sum of all cross products around loop adds area for
1356  * edges away from origin and subtracts area for edges
1357  * near origin, leaving twice the area of the polygon.
1358  */
1359  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
1360  struct edgeuse *next_eu;
1361  struct vertex *vp, *next_vp;
1362  vect_t cross;
1363  point_t p1, p2;
1364 
1365  vp = eu->vu_p->v_p;
1366 
1367  next_eu = BU_LIST_PNEXT_CIRC(edgeuse, &eu->l);
1368  NMG_CK_EDGEUSE(next_eu);
1369  NMG_CK_VERTEXUSE(next_eu->vu_p);
1370 
1371  next_vp = next_eu->vu_p->v_p;
1372  NMG_CK_VERTEX(next_vp);
1373  NMG_CK_VERTEX_G(next_vp->vg_p);
1374 
1375  VSUB2(p1, vp->vg_p->coord, trans);
1376  VSUB2(p2, next_vp->vg_p->coord, trans);
1377  VCROSS(cross, p1, p2);
1378  VADD2(plane, plane, cross);
1379  }
1380 
1381  area = 0.5 * MAGNITUDE(plane);
1382 
1383  /* Error if the area is too small to unitize the normal vector */
1384  if (MAGNITUDE(plane) < VDIVIDE_TOL) return -1.0;
1385  VUNITIZE(plane);
1386 
1387  /* calculate plane[W] as average distance to plane */
1388  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
1389  pt_dot_plane += VDOT(plane, eu->vu_p->v_p->vg_p->coord);
1390  pt_count++;
1391  }
1392 
1393  /* Error if we don't have at least 3 vertices in this loop */
1394  if (pt_count < 3)
1395  return (fastf_t)(-1.0);
1396 
1397  plane[W] = pt_dot_plane/pt_count;
1398  HMOVE(pl, plane);
1399 
1400  return area;
1401 }
1402 
1403 
1404 /**
1405  * Calculates a plane equation and the area of a loop
1406  *
1407  * This function only works correctly when the loop represents a
1408  * "simple polygon" (i.e. simple meaning not self intersecting) and
1409  * the loop is "not" a "weakly simple polygon" (i.e. weakly simple
1410  * meaning is simple and at least one vertex is reused by the polygon).
1411  *
1412  * returns:
1413  * the area of the loop
1414  * positive area is ccw loop
1415  * negative area is cw loop
1416  * zero is no area or degenerate loop
1417  *
1418  * pl is assigned the plane equation for the loop
1419  *
1420  * NOTE: The rotation cw/ccw is actual, meaning not relative to the
1421  * faceuse normal.
1422  */
1423 fastf_t
1424 nmg_loop_plane_area2(const struct loopuse *lu, fastf_t *pl, const struct bn_tol *tol)
1425 {
1426  struct edgeuse *eu;
1427  size_t cnt;
1428  vect_t zaxis = {0.0, 0.0, 1.0};
1429  vect_t sum, cog;
1430  mat_t mat;
1431  point_t pt, pt_next, pt_1st;
1432  fastf_t scale;
1433  fastf_t area = 0.0;
1434  fastf_t min, abs_coord, mag;
1435  fastf_t *tmp;
1436  int i;
1437 
1438  NMG_CK_LOOPUSE(lu);
1439 
1440  HSETALL(pl, 0.0);
1441 
1442  if (UNLIKELY(BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)) {
1443  area = 0.0;
1444  goto out;
1445  }
1446 
1447  cnt = 0;
1448  min = MAX_FASTF;
1449  /* find the smallest floating point value which will be used
1450  * to determine the scaling factor, which is done to improve
1451  * floating point precision
1452  */
1453  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
1454  tmp = eu->vu_p->v_p->vg_p->coord;
1455  for (i = 0; i < 3 ; i++) {
1456  abs_coord = fabs(tmp[i]);
1457  if ((abs_coord) < min && (abs_coord > SMALL_FASTF)) {
1458  min = abs_coord;
1459  }
1460  }
1461  cnt++;
1462  }
1463 
1464  if (cnt < 3) {
1465  area = 0.0;
1466  goto out;
1467  }
1468 
1469  scale = 0.5 / min;
1470  if (scale < 1.0) {
1471  scale = 1.0;
1472  } else if (scale > (1.0 / tol->dist)) {
1473  scale = 1.0 / tol->dist;
1474  }
1475 
1476  /* find and scale first vertex of polygon */
1477  VMOVE(pt_1st, (BU_LIST_FIRST(edgeuse, &eu->l))->vu_p->v_p->vg_p->coord);
1478  VINTCLAMP(pt_1st);
1479  VSCALE(pt_1st, pt_1st, scale);
1480 
1481  VSETALL(sum, 0.0);
1482  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
1483  NMG_CK_EDGEUSE(eu);
1484  /* scale up the coordinates to reduce floating point error */
1485  VSCALE(pt, eu->vu_p->v_p->vg_p->coord, scale);
1486  VSCALE(pt_next, (BU_LIST_PNEXT_CIRC(edgeuse, &eu->l))->vu_p->v_p->vg_p->coord, scale);
1487  /* move polygon to origin */
1488  VSUB2(pt, pt, pt_1st);
1489  VSUB2(pt_next, pt_next, pt_1st);
1490  VINTCLAMP(pt);
1491  VINTCLAMP(pt_next);
1492  pl[X] += (EQUAL(pt[Y], pt_next[Y]) ? (0.0) : (pt[Y] - pt_next[Y])) * (pt[Z] + pt_next[Z]);
1493  pl[Y] += (EQUAL(pt[Z], pt_next[Z]) ? (0.0) : (pt[Z] - pt_next[Z])) * (pt[X] + pt_next[X]);
1494  pl[Z] += (EQUAL(pt[X], pt_next[X]) ? (0.0) : (pt[X] - pt_next[X])) * (pt[Y] + pt_next[Y]);
1495  VADD2(sum, sum, pt);
1496  }
1497 
1498  if (ZERO(pl[X]) && ZERO(pl[Y]) && ZERO(pl[Z])) {
1499  HSETALL(pl, 0.0);
1500  area = 0.0;
1501  goto out;
1502  }
1503 
1504  VINTCLAMP(pl);
1505 
1506  /* undo scaling of center of polygon (i.e. center of gravity) */
1507  VSCALE(cog, sum, 1.0/(cnt * scale));
1508 
1509  /* Unitize plane. If magnitude is zero, the plane is undefined and
1510  * we can not compute the area with this function.
1511  */
1512  mag = MAGNITUDE(pl);
1513  if (ZERO(mag)) {
1514  HSETALL(pl, 0.0);
1515  area = 0.0;
1516  goto out;
1517  }
1518  VSCALE(pl, pl, 1.0/mag);
1519  pl[H] = VDOT(cog, pl);
1520 
1521  /* find rotation matrix of polygon normal to z-axis */
1522  bn_mat_fromto(mat, (const fastf_t *)pl, (const fastf_t *)zaxis, tol);
1523 
1524  /* compute area of polygon */
1525  area = 0.0;
1526  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
1527  /* rotate coordinates to xy plane */
1528  MAT4X3PNT(pt, mat, eu->vu_p->v_p->vg_p->coord);
1529  MAT4X3PNT(pt_next, mat, (BU_LIST_PNEXT_CIRC(edgeuse, &eu->l))->vu_p->v_p->vg_p->coord);
1530  area += ((pt[X]*pt_next[Y]) - (pt_next[X]*pt[Y]));
1531  }
1532  area = area / 2.0;
1533 
1534 out:
1535  return area;
1536 }
1537 
1538 
1539 /**
1540  * Calculate face geometry using a least squares fit or Newell's
1541  * method.
1542  *
1543  * If fu does not already have a face_g_plane associated, only
1544  * vertices in fu will participate, and if it has only one loop
1545  * Newell's method will be used rather than a least square fit.
1546  *
1547  * if fu has a face_g_plane, then all vertices in any face that
1548  * references the same face_g_plane will participate in the fit for
1549  * the face plane.
1550  *
1551  * Returns:
1552  * 0 - All is well
1553  * 1 - Failed to calculate face geometry
1554  *
1555  */
1556 int
1557 nmg_calc_face_plane(struct faceuse *fu_in, fastf_t *pl)
1558 {
1559  double one_over_vertex_count;
1560  fastf_t det;
1561  fastf_t max_dist=(-MAX_FASTF);
1562  fastf_t min_dist=MAX_FASTF;
1563  int failed=0;
1564  int got_dir=0;
1565  int i;
1566  int loop_count=0;
1567  mat_t inverse = MAT_INIT_ZERO;
1568  mat_t matrix = MAT_INIT_ZERO;
1569  plane_t old_pl = HINIT_ZERO;
1570  struct bu_ptbl verts;
1571  struct face *f;
1572  struct face_g_plane *fg;
1573  struct faceuse *fu;
1574  struct loopuse *lu;
1575  vect_t vsum = VINIT_ZERO;
1576 
1577  fu = fu_in;
1578  NMG_CK_FACEUSE(fu);
1579 
1580  HSETALL(pl, 0.0); /* sanity */
1581 
1582  /* find an OT_SAME loop to use for calculating general direction of normal */
1583  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
1584  if (!got_dir && BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_EDGEUSE_MAGIC) {
1585  /* get general direction for face normal */
1586  nmg_loop_plane_newell(lu, old_pl);
1587  if (!ZERO(old_pl[X]) || !ZERO(old_pl[Y]) || !ZERO(old_pl[Z]))
1588  got_dir = 1;
1589  }
1590 
1591  loop_count++;
1592  }
1593 
1594  if (!got_dir) {
1595  failed = 1;
1596  goto out;
1597  }
1598 
1599  f = fu->f_p;
1600  NMG_CK_FACE(f);
1601  fg = f->g.plane_p;
1602  if (fg) {
1603  struct face *f1;
1604 
1605  NMG_CK_FACE_G_PLANE(fg);
1606 
1607  /* count loops using this face geometry */
1608  loop_count = 0;
1609  for (BU_LIST_FOR (f1, face, &fg->f_hd)) {
1610  for (BU_LIST_FOR (lu, loopuse, &f1->fu_p->lu_hd))
1611  loop_count++;
1612  }
1613 
1614  /* if this face geometry only has one loop using it, just use Newell's method */
1615  if (loop_count < 2) {
1616  HMOVE(pl, old_pl);
1617  failed = 0;
1618  goto out;
1619  }
1620 
1621  nmg_tabulate_face_g_verts(&verts, fg);
1622  } else {
1623  /* just use vertices from this faceuse */
1624  /* If this faceuse only has one loop, just use Newell's method */
1625  if (loop_count < 2) {
1626  HMOVE(pl, old_pl);
1627  failed = 0;
1628  goto out;
1629  }
1630 
1631  nmg_vertex_tabulate(&verts, &fu->l.magic);
1632  }
1633 
1634  /* Get the direction for the plane normal in "old_pl".
1635  * Make sure we are dealing with OT_SAME or OT_UNSPEC faceuse.
1636  */
1637  if (fu->orientation != OT_UNSPEC) {
1638  if (fu->orientation != OT_SAME)
1639  fu = fu->fumate_p;
1640  if (fu->orientation != OT_SAME) {
1641  bu_log("nmg_calc_face_plane: fu %p has no OT_SAME use\n", (void *)fu);
1642  bu_ptbl_free(&verts);
1643  failed = 1;
1644  goto out;
1645  }
1646  }
1647 
1648  /* build matrix */
1649  one_over_vertex_count = 1.0/(double)(BU_PTBL_END(&verts));
1650 
1651  for (i = 0; i < BU_PTBL_END(&verts); i++) {
1652  struct vertex *v;
1653  struct vertex_g *vg;
1654 
1655  v = (struct vertex *)BU_PTBL_GET(&verts, i);
1656  vg = v->vg_p;
1657 
1658  matrix[0] += vg->coord[X] * vg->coord[X];
1659  matrix[1] += vg->coord[X] * vg->coord[Y];
1660  matrix[2] += vg->coord[X] * vg->coord[Z];
1661  matrix[5] += vg->coord[Y] * vg->coord[Y];
1662  matrix[6] += vg->coord[Y] * vg->coord[Z];
1663  matrix[10] += vg->coord[Z] * vg->coord[Z];
1664 
1665  vsum[X] += vg->coord[X];
1666  vsum[Y] += vg->coord[Y];
1667  vsum[Z] += vg->coord[Z];
1668  }
1669  matrix[4] = matrix[1];
1670  matrix[8] = matrix[2];
1671  matrix[9] = matrix[6];
1672  matrix[15] = 1.0;
1673 
1674 
1675  /* Check that we don't have a singular matrix */
1676  det = bn_mat_determinant(matrix);
1677 
1678  if (!ZERO(det)) {
1679  fastf_t inv_len_pl;
1680 
1681  /* invert matrix */
1682  bn_mat_inv(inverse, matrix);
1683 
1684  /* get normal vector */
1685  MAT4X3PNT(pl, inverse, vsum);
1686 
1687  /* unitize direction vector */
1688  inv_len_pl = 1.0/(MAGNITUDE(pl));
1689  HSCALE(pl, pl, inv_len_pl);
1690 
1691  /* get average vertex coordinates */
1692  VSCALE(vsum, vsum, one_over_vertex_count);
1693 
1694  /* get distance from plane to origin */
1695  pl[H] = VDOT(pl, vsum);
1696 
1697  /* make sure it points in the correct direction */
1698  if (VDOT(pl, old_pl) < -SMALL_FASTF) {
1699  HREVERSE(pl, pl);
1700  }
1701  } else {
1702  struct vertex *v, *v0;
1703  int x_same=1;
1704  int y_same=1;
1705  int z_same=1;
1706 
1707  /* singular matrix, may occur if all vertices have the same zero
1708  * component.
1709  */
1710  v0 = (struct vertex *)BU_PTBL_GET(&verts, 0);
1711  for (i=1; i<BU_PTBL_END(&verts); i++) {
1712  v = (struct vertex *)BU_PTBL_GET(&verts, i);
1713 
1714  if (!ZERO(v->vg_p->coord[X] - v0->vg_p->coord[X]))
1715  x_same = 0;
1716  if (!ZERO(v->vg_p->coord[Y] - v0->vg_p->coord[Y]))
1717  y_same = 0;
1718  if (!ZERO(v->vg_p->coord[Z] - v0->vg_p->coord[Z]))
1719  z_same = 0;
1720 
1721  if (!x_same && !y_same && !z_same)
1722  break;
1723  }
1724 
1725  if (x_same) {
1726  HSET(pl, 1.0, 0.0, 0.0, 0.0);
1727  } else if (y_same) {
1728  HSET(pl, 0.0, 1.0, 0.0, 0.0);
1729  } else if (z_same) {
1730  HSET(pl, 0.0, 0.0, 1.0, 0.0);
1731  }
1732 
1733  if (x_same || y_same || z_same) {
1734  /* get average vertex coordinates */
1735  VSCALE(vsum, vsum, one_over_vertex_count);
1736 
1737  /* get distance from plane to origin */
1738  pl[H] = VDOT(pl, vsum);
1739 
1740  /* make sure it points in the correct direction */
1741  if (VDOT(pl, old_pl) < -SMALL_FASTF) {
1742  HREVERSE(pl, pl);
1743  }
1744  } else {
1745  bu_log("nmg_calc_face_plane: Cannot calculate plane for fu %p\n", (void *)fu);
1746  nmg_pr_fu_briefly(fu, (char *)NULL);
1747  bu_log("%ld verts\n", BU_PTBL_END(&verts));
1748  failed = 1;
1749  goto out;
1750  }
1751  }
1752 
1753  /* make sure plane is at center of range of vertices */
1754  for (i=0; i<BU_PTBL_END(&verts); i++) {
1755  struct vertex *v;
1756  struct vertex_g *vg;
1757  fastf_t dist;
1758 
1759  v = (struct vertex *)BU_PTBL_GET(&verts, i);
1760  vg = v->vg_p;
1761 
1762  dist = DIST_PT_PLANE(vg->coord, pl);
1763  if (dist > max_dist)
1764  max_dist = dist;
1765  if (dist < min_dist)
1766  min_dist = dist;
1767  }
1768 
1769  pl[H] += (max_dist + min_dist)/2.0;
1770 
1771  bu_ptbl_free(&verts);
1772 
1773  failed = 0;
1774 
1775 out:
1776  return failed;
1777 
1778 }
1779 
1780 
1781 /**
1782  * interface to nmg_calc_face_plane(), calls nmg_face_g with the
1783  * resulting plane
1784  */
1785 int
1786 nmg_calc_face_g(struct faceuse *fu)
1787 {
1788  plane_t pl;
1789  int ret_val;
1790 
1791  ret_val = nmg_calc_face_plane(fu, pl);
1792 
1793  if (!ret_val)
1794  nmg_face_g(fu, pl);
1795 
1796  return ret_val;
1797 }
1798 
1799 
1800 /**
1801  * The following routines calculate surface area of NMG objects. Note
1802  * that this includes all surfaces, not just external surfaces, i.e.,
1803  * an NMG object consisting of two adjacent cubes with a coincident
1804  * face will have a surface area of 12*s*s (s is length of one side)
1805  */
1806 fastf_t
1807 nmg_faceuse_area(const struct faceuse *fu)
1808 {
1809  struct loopuse *lu;
1810  plane_t plane;
1811  fastf_t area=0.0;
1812  fastf_t tmp_area;
1813 
1814  NMG_CK_FACEUSE(fu);
1815 
1816  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
1817  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
1818  continue;
1819 
1820  tmp_area = nmg_loop_plane_area(lu, plane);
1821  if (tmp_area < 0.0)
1822  continue;
1823 
1824  if (lu->orientation == OT_SAME)
1825  area += tmp_area;
1826  else if (lu->orientation == OT_OPPOSITE)
1827  area -= tmp_area;
1828  else {
1829  bu_log("nmg_faceuse_area: Cannot calculate area (lu with %s orientation)\n",
1830  nmg_orientation(lu->orientation));
1831  return (fastf_t)-1.0;
1832  }
1833  }
1834 
1835  return area;
1836 }
1837 
1838 
1839 fastf_t
1840 nmg_shell_area(const struct shell *s)
1841 {
1842  fastf_t area=0.0;
1843  fastf_t tmp_area;
1844  struct faceuse *fu;
1845 
1846  NMG_CK_SHELL(s);
1847 
1848  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
1849  if (fu->orientation != OT_SAME)
1850  continue;
1851 
1852  tmp_area = nmg_faceuse_area(fu);
1853  if (tmp_area < 0.0)
1854  continue;
1855 
1856  area += tmp_area;
1857  }
1858 
1859  return area;
1860 }
1861 
1862 
1863 fastf_t
1864 nmg_region_area(const struct nmgregion *r)
1865 {
1866  struct shell *s;
1867  fastf_t area=0.0;
1868 
1869  NMG_CK_REGION(r);
1870 
1871  for (BU_LIST_FOR (s, shell, &r->s_hd))
1872  area += nmg_shell_area(s);
1873 
1874  return area;
1875 }
1876 
1877 
1878 fastf_t
1879 nmg_model_area(const struct model *m)
1880 {
1881  struct nmgregion *r;
1882  fastf_t area=0.0;
1883 
1884  NMG_CK_MODEL(m);
1885 
1886  for (BU_LIST_FOR (r, nmgregion, &m->r_hd))
1887  area += nmg_region_area(r);
1888 
1889  return area;
1890 }
1891 
1892 
1893 /**
1894  * Make sure that the list of intersection points doesn't contain any
1895  * vertexuses from loops whose bounding boxes don;t overlap the
1896  * bounding box of a loop in the given faceuse.
1897  *
1898  * This is really a special purpose routine to help the intersection
1899  * operations of the boolean process. The only reason it's here
1900  * instead of nmg_inter.c is that it knows too much about the format
1901  * and contents of an bu_ptbl structure.
1902  */
1903 void
1904 nmg_purge_unwanted_intersection_points(struct bu_ptbl *vert_list, fastf_t *mag_list, const struct faceuse *fu, const struct bn_tol *tol)
1905 {
1906  int i;
1907  int j;
1908  struct vertexuse *vu;
1909  struct loopuse *lu;
1910  const struct loop_g *lg;
1911  const struct loopuse *fu2lu;
1912  const struct loop_g *fu2lg = (const struct loop_g *)NULL;
1913  int overlap = 0;
1914 
1915  NMG_CK_FACEUSE(fu);
1916  BN_CK_TOL(tol);
1917 
1918  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1919  bu_log("nmg_purge_unwanted_intersection_points(%p, %p)\n", (void *)vert_list, (void *)fu);
1920  bu_log("\t%ld vertexuses in list\n", vert_list->end);
1921  }
1922 
1923  for (i = 0; i < vert_list->end; i++) {
1924  vu = (struct vertexuse *)vert_list->buffer[i];
1925  if (vu->l.magic != NMG_VERTEXUSE_MAGIC) {
1926  /* vertexuse probably killed by nmg_repair_v_near_v() */
1927  /* delete the entry from the vertex list */
1928  for (j=i; j < vert_list->end; j++) {
1929  vert_list->buffer[j] = vert_list->buffer[j+1];
1930  mag_list[j] = mag_list[j+1];
1931  }
1932 
1933  --(vert_list->end);
1934  vert_list->buffer[vert_list->end] = (long *)NULL;
1935  mag_list[vert_list->end] = MAX_FASTF;
1936  --i;
1937  continue;
1938  }
1939  NMG_CK_VERTEXUSE(vu);
1940  lu = nmg_find_lu_of_vu(vu);
1941  NMG_CK_LOOPUSE(lu);
1942  lg = lu->l_p->lg_p;
1943  NMG_CK_LOOP_G(lg);
1944 
1945  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1946  bu_log("vu[%d]: %p (%g %g %g) lu: %p %s\n",
1947  i,
1948  (void *)vu, V3ARGS(vu->v_p->vg_p->coord),
1949  (void *)lu, nmg_orientation(lu->orientation));
1950  bu_log("\tlu BBox: (%g %g %g) (%g %g %g)\n",
1951  V3ARGS(lg->min_pt), V3ARGS(lg->max_pt));
1952  }
1953  if (lu->up.fu_p->f_p == fu->f_p)
1954  bu_log("I'm checking against my own face?\n");
1955 
1956  /* If the bounding box of a loop doesn't intersect the
1957  * bounding box of a loop in the other face, it shouldn't
1958  * be on the list of intersecting loops.
1959  */
1960  overlap = 0;
1961  for (BU_LIST_FOR (fu2lu, loopuse, &fu->lu_hd)) {
1962  NMG_CK_LOOPUSE(fu2lu);
1963  NMG_CK_LOOP(fu2lu->l_p);
1964 
1965  switch (fu2lu->orientation) {
1966  case OT_BOOLPLACE:
1967  /* If this loop is destined for removal
1968  * by sanitize(), skip it.
1969  */
1970  continue;
1971  case OT_UNSPEC:
1972  /* If this is a loop of a lone vertex,
1973  * it was deposited into the
1974  * other loop as part of an intersection
1975  * operation, and is quite important.
1976  */
1977  if (BU_LIST_FIRST_MAGIC(&fu2lu->down_hd) != NMG_VERTEXUSE_MAGIC)
1978  bu_log("nmg_purge_unwanted_intersection_points() non self-loop OT_UNSPEC vertexuse in fu2?\n");
1979  break;
1980  case OT_SAME:
1981  case OT_OPPOSITE:
1982  break;
1983  default:
1984  bu_log("nmg_purge_unwanted_intersection_points encountered %s loop in fu2\n",
1985  nmg_orientation(fu2lu->orientation));
1986  /* Process it anyway */
1987  break;
1988  }
1989 
1990  fu2lg = fu2lu->l_p->lg_p;
1991  NMG_CK_LOOP_G(fu2lg);
1992 
1993  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1994  bu_log("\tfu2lu BBox: (%g %g %g) (%g %g %g) %s\n",
1995  V3ARGS(fu2lg->min_pt), V3ARGS(fu2lg->max_pt),
1996  nmg_orientation(fu2lu->orientation));
1997  }
1998 
1999  if (V3RPP_OVERLAP_TOL(fu2lg->min_pt, fu2lg->max_pt, lg->min_pt, lg->max_pt, tol->dist)) {
2000  overlap = 1;
2001  break;
2002  }
2003  }
2004  if (!overlap) {
2005  /* why is this vertexuse in the list? */
2006  if (RTG.NMG_debug & DEBUG_POLYSECT) {
2007  bu_log("nmg_purge_unwanted_intersection_points This little bugger slipped in somehow. Deleting it from the list.\n");
2008  nmg_pr_vu_briefly(vu, (char *)NULL);
2009  }
2010  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC &&
2011  lu->orientation == OT_UNSPEC) {
2012  /* Drop this loop of a single vertex in sanitize() */
2013  if (RTG.NMG_debug & DEBUG_POLYSECT)
2014  bu_log("nmg_purge_unwanted_intersection_points() remarking as OT_BOOLPLACE\n");
2015  lu->orientation =
2016  lu->lumate_p->orientation = OT_BOOLPLACE;
2017  }
2018 
2019  /* delete the entry from the vertex list */
2020  for (j=i; j < vert_list->end; j++) {
2021  vert_list->buffer[j] = vert_list->buffer[j+1];
2022  mag_list[j] = mag_list[j+1];
2023  }
2024 
2025  --(vert_list->end);
2026  vert_list->buffer[vert_list->end] = (long *)NULL;
2027  mag_list[vert_list->end] = MAX_FASTF;
2028  --i;
2029  }
2030  }
2031  if (RTG.NMG_debug & DEBUG_POLYSECT)
2032  bu_log("\tAt end of nmg_purge_unwanted_intersection_points, %ld vertexuses in list\n",
2033  vert_list->end);
2034 }
2035 
2036 
2037 /**
2038  * if the given vertexuse "vu" is in the table given by "b" or if "vu"
2039  * references a vertex which is referenced by a vertexuse in the table,
2040  * then we return 1. Otherwise, we return 0.
2041  */
2042 int
2043 nmg_in_or_ref(struct vertexuse *vu, struct bu_ptbl *b)
2044 {
2045  union {
2046  struct vertexuse **vu;
2047  uint32_t **magic_p;
2048  } p;
2049  int i;
2050 
2051  p.magic_p = (uint32_t **)b->buffer;
2052  for (i=0; i < b->end; ++i) {
2053  if (p.vu[i] && *p.magic_p[i] == NMG_VERTEXUSE_MAGIC &&
2054  (p.vu[i] == vu || p.vu[i]->v_p == vu->v_p))
2055  return 1;
2056  }
2057  return 0;
2058 }
2059 
2060 
2061 /**
2062  * Re-compute all the bounding boxes in the NMG model.
2063  * Bounding boxes are presently found in these structures:
2064  * loop_g
2065  * face_g
2066  * shell_a
2067  * nmgregion_a
2068  * The re-bounding must be performed in a bottom-up manner,
2069  * computing the loops first, and working up to the nmgregions.
2070  */
2071 void
2072 nmg_rebound(struct model *m, const struct bn_tol *tol)
2073 {
2074  struct nmgregion *r;
2075  struct shell *s;
2076  struct faceuse *fu;
2077  struct face *f;
2078  struct loopuse *lu;
2079  struct loop *l;
2080  register int *flags;
2081 
2082  NMG_CK_MODEL(m);
2083  BN_CK_TOL(tol);
2084 
2085  flags = (int *)bu_calloc(m->maxindex*2, sizeof(int), "rebound flags[]");
2086 
2087  for (BU_LIST_FOR (r, nmgregion, &m->r_hd)) {
2088  NMG_CK_REGION(r);
2089  for (BU_LIST_FOR (s, shell, &r->s_hd)) {
2090  NMG_CK_SHELL(s);
2091 
2092  /* Loops in faces in shell */
2093  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
2094  NMG_CK_FACEUSE(fu);
2095  /* Loops in face */
2096  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
2097  NMG_CK_LOOPUSE(lu);
2098  l = lu->l_p;
2099  NMG_CK_LOOP(l);
2100  if (NMG_INDEX_FIRST_TIME(flags, l))
2101  nmg_loop_g(l, tol);
2102  }
2103  }
2104  /* Faces in shell */
2105  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
2106  NMG_CK_FACEUSE(fu);
2107  f = fu->f_p;
2108  NMG_CK_FACE(f);
2109 
2110  /* Rebound the face */
2111  if (NMG_INDEX_FIRST_TIME(flags, f))
2112  nmg_face_bb(f, tol);
2113  }
2114 
2115  /* Wire loops in shell */
2116  for (BU_LIST_FOR (lu, loopuse, &s->lu_hd)) {
2117  NMG_CK_LOOPUSE(lu);
2118  l = lu->l_p;
2119  NMG_CK_LOOP(l);
2120  if (NMG_INDEX_FIRST_TIME(flags, l))
2121  nmg_loop_g(l, tol);
2122  }
2123 
2124  /*
2125  * Rebound the shell.
2126  * This routine handles wire edges and lone vertices.
2127  */
2128  if (NMG_INDEX_FIRST_TIME(flags, s))
2129  nmg_shell_a(s, tol);
2130  }
2131 
2132  /* Rebound the nmgregion */
2133  nmg_region_a(r, tol);
2134  }
2135 
2136  bu_free((char *)flags, "rebound flags[]");
2137 }
2138 
2139 
2140 void
2141 nmg_count_shell_kids(const struct model *m, size_t *total_faces, size_t *total_wires, size_t *total_points)
2142 {
2143  short *tbl;
2144 
2145  const struct nmgregion *r;
2146  const struct shell *s;
2147  const struct faceuse *fu;
2148  const struct loopuse *lu;
2149  const struct edgeuse *eu;
2150 
2151  NMG_CK_MODEL(m);
2152 
2153  tbl = (short *)bu_calloc(m->maxindex+1, sizeof(char),
2154  "face/wire/point counted table");
2155 
2156  *total_faces = *total_wires = *total_points = 0;
2157  for (BU_LIST_FOR (r, nmgregion, &m->r_hd)) {
2158  for (BU_LIST_FOR (s, shell, &r->s_hd)) {
2159  if (s->vu_p) {
2160  (*total_points)++;
2161  continue;
2162  }
2163  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
2164  if (NMG_INDEX_TEST_AND_SET(tbl, fu->f_p))
2165  (*total_faces)++;
2166  }
2167  for (BU_LIST_FOR (lu, loopuse, &s->lu_hd)) {
2168  if (NMG_INDEX_TEST_AND_SET(tbl, lu->l_p))
2169  (*total_wires)++;
2170  }
2171  for (BU_LIST_FOR (eu, edgeuse, &s->eu_hd)) {
2172  if (NMG_INDEX_TEST_AND_SET(tbl, eu->e_p))
2173  (*total_wires)++;
2174  }
2175  }
2176  }
2177 
2178  bu_free((char *)tbl, "face/wire/point counted table");
2179 }
2180 
2181 
2182 /**
2183  * private support routine for nmg_close_shell creates an array of
2184  * indices into a table of edgeuses, ordered end-to-end. This may or
2185  * may not create an actual loop.
2186  *
2187  * Arguments:
2188  * tbl is the table (provided by caller)
2189  * index is the array of indices created by order_tbl
2190  * tbl_size is the size of the table (provided by caller)
2191  * loop_size is the number of edgeuses in the loop (calculated by order_tbl)
2192  */
2193 HIDDEN void
2194 order_tbl(struct bu_ptbl *tbl, int start_idx, int **idx, int tbl_size, int *loop_size)
2195 {
2196  int i, j, k;
2197  int found;
2198  struct edgeuse *eu, *eu1;
2199  struct vertex *start_v;
2200 
2201  /* create an index into the table, ordered to create a loop */
2202  if (*idx == NULL)
2203  (*idx) = (int *)bu_calloc(tbl_size, sizeof(int), "Table index");
2204 
2205  for (i=0; i<tbl_size; i++)
2206  (*idx)[i] = (-1);
2207 
2208  /* start the loop at idx = start_idx */
2209  (*idx)[0] = start_idx;
2210  *loop_size = 1;
2211  eu = (struct edgeuse *)BU_PTBL_GET(tbl, start_idx);
2212  start_v = eu->vu_p->v_p;
2213  i = 0;
2214  found = 1;
2215  while (eu->eumate_p->vu_p->v_p != start_v && found) {
2216  found = 0;
2217 
2218  /* Look for edgeuse that starts where "eu" ends */
2219  for (j=0; j<tbl_size; j++) {
2220  int already_used = 0;
2221 
2222  eu1 = (struct edgeuse *)BU_PTBL_GET(tbl, j);
2223 
2224  /* don't use same edgeuse twice! */
2225  for (k=0; k<(*loop_size); k++) {
2226  if (eu1 == (struct edgeuse *)BU_PTBL_GET(tbl, (*idx)[k])) {
2227  already_used = 1;
2228  break;
2229  }
2230  }
2231  if (already_used)
2232  continue;
2233  if (eu1->vu_p->v_p == eu->eumate_p->vu_p->v_p) {
2234  /* Found it */
2235  found = 1;
2236  (*idx)[++i] = j;
2237  (*loop_size)++;
2238  eu = eu1;
2239  break;
2240  }
2241  }
2242  }
2243 }
2244 
2245 
2246 /**
2247  * Examines the passed shell and, if there are holes, closes them
2248  * note that not much care is taken as to how the holes are closed
2249  * so the results are not entirely predictable.
2250  * A list of free edges is created (edges bounding only one face).
2251  * New faces are constructed by taking two consecutive edges
2252  * and making a face. The newly created edge is added to the list
2253  * of free edges and the two used ones are removed.
2254  */
2255 void
2256 nmg_close_shell(struct shell *s, const struct bn_tol *tol)
2257 {
2258  struct bu_ptbl eu_tbl; /* table of free edgeuses from shell */
2259  struct bu_ptbl vert_tbl; /* table of vertices for use in nmg_cface */
2260  int *idx; /* array of indices into eu_tbl, ordered to form a loop */
2261  int loop_size; /* number of edgeuses in loop */
2262  struct faceuse *fu;
2263  struct loopuse *lu;
2264  struct edgeuse *eu;
2265  int start_loop;
2266  int i;
2267  int found;
2268 
2269  if (RTG.NMG_debug & DEBUG_BASIC)
2270  bu_log("nmg_close_shell: s = %p\n", (void *)s);
2271 
2272  NMG_CK_SHELL(s);
2273  BN_CK_TOL(tol);
2274 
2275  idx = NULL;
2276 
2277  /* construct a table of free edges */
2278  (void)bu_ptbl_init(&eu_tbl, 64, " &eu_tbl ");
2279 
2280  /* loop through all the faces in the shell */
2281  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
2282  NMG_CK_FACEUSE(fu);
2283  /* only look at OT_SAME faces */
2284  if (fu->orientation == OT_SAME) {
2285  /* loop through each loopuse in the face */
2286  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
2287  NMG_CK_LOOPUSE(lu);
2288  /* ignore loops that are just a vertex */
2289  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) ==
2291  continue;
2292 
2293  /* loop through all the edgeuses in the loop */
2294  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
2295  NMG_CK_EDGEUSE(eu);
2296  /* if this edgeuse is a free edge, add its mate to the list */
2297  if (eu->radial_p == eu->eumate_p)
2298  (void)bu_ptbl_ins(&eu_tbl, (long *) eu->eumate_p);
2299  }
2300  }
2301  }
2302  }
2303 
2304  /* if there is nothing in our list of free edges, the shell is already closed */
2305  if (BU_PTBL_END(&eu_tbl) == 0) {
2306  bu_ptbl_free(&eu_tbl);
2307  return;
2308  }
2309 
2310  /* create a table of vertices */
2311  (void)bu_ptbl_init(&vert_tbl, 64, " &vert_tbl ");
2312 
2313  while (BU_PTBL_END(&eu_tbl)) {
2314  int give_up_on_face=0;
2315 
2316  /* Create an index into the table that orders the edgeuses into a loop */
2317  start_loop = -1;
2318  loop_size = 0;
2319  while (loop_size < 2) {
2320  if (++start_loop > BU_PTBL_END(&eu_tbl) - 3)
2321  break;
2322  order_tbl(&eu_tbl, start_loop, &idx, BU_PTBL_END(&eu_tbl), &loop_size);
2323  }
2324 
2325  /* Create new faces to close the shell */
2326  while (loop_size > 3) {
2327  struct edgeuse *eu1, *eu2=NULL, *eu_new = NULL;
2328  struct edgeuse **eu_used; /* array of edgeuses used, for deletion */
2329  int edges_used=0; /* number of edges used in making a face */
2330  int found_face=0; /* flag - indicates that a face with the correct normal will be created */
2331  int start_idx, end_idx; /* start and stop index for loop */
2332  int coplanar; /* flag - indicates entire loop is coplanar */
2333  plane_t pl1, pl2; /* planes for checking coplanarity of loop */
2334  point_t pt[3]; /* points for calculating planes */
2335 
2336  /* Look for an easy way out, maybe this loop is planar */
2337  /* first, calculate a plane from the first three non-collinear vertices */
2338  start_idx = 0;
2339  end_idx = start_idx + 3;
2340 
2341  for (i=start_idx; i<end_idx; i++) {
2342  eu = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[i]);
2343  VMOVE(pt[i-start_idx], eu->vu_p->v_p->vg_p->coord);
2344  }
2345  while (bn_mk_plane_3pts(pl1, pt[0], pt[1], pt[2], tol) && end_idx<loop_size) {
2346  start_idx++;
2347  end_idx++;
2348  for (i=start_idx; i<end_idx; i++) {
2349  eu = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[i]);
2350  VMOVE(pt[i-start_idx], eu->vu_p->v_p->vg_p->coord);
2351  }
2352  }
2353  if (end_idx == loop_size) {
2354  if (BU_PTBL_END(&eu_tbl) > loop_size) {
2355  int old_loop_size=loop_size;
2356 
2357  loop_size = 0;
2358  while (loop_size < 3) {
2359  found = 1;
2360  for (start_loop=idx[0]+1; start_loop < BU_PTBL_END(&eu_tbl)-3; start_loop++) {
2361  found = 0;
2362  for (i=0; i<old_loop_size; i++) {
2363  if (idx[i] == start_loop) {
2364  found = 1;
2365  break;
2366  }
2367  }
2368  if (!found)
2369  break;
2370  }
2371  if (found) {
2372  /* Could not even make a plane, this is some serious screw-up */
2373  bu_log("nmg_close_shell: cannot make any planes from loop, old_loop_size = %d\n", old_loop_size);
2374  for (i=0; i<old_loop_size; i++) {
2375  eu= (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[i]);
2376  bu_log("(%g %g %g) <-> (%g %g %g)\n",
2377  V3ARGS(eu->vu_p->v_p->vg_p->coord),
2378  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord));
2379  }
2380  bu_bomb("nmg_close_shell: cannot make any planes from loop\n");
2381  }
2382  order_tbl(&eu_tbl, start_loop, &idx, BU_PTBL_END(&eu_tbl), &loop_size);
2383  }
2384  }
2385  }
2386 
2387  /* now we have one plane, let's check the others */
2388  coplanar = 1;
2389  while (end_idx < loop_size && coplanar) {
2390  end_idx +=3;
2391  if (end_idx > loop_size)
2392  end_idx = loop_size;
2393  start_idx = end_idx - 3;
2394 
2395  for (i=start_idx; i<end_idx; i++) {
2396  eu = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[i]);
2397  VMOVE(pt[i-start_idx], eu->vu_p->v_p->vg_p->coord);
2398  }
2399 
2400  /* if these three points make a plane, is it coplanar with
2401  * our first one??? */
2402  if (!bn_mk_plane_3pts(pl2, pt[0], pt[1], pt[2], tol)) {
2403  if ((i=bn_coplanar(pl1, pl2, tol)) < 1)
2404  coplanar = 0;
2405  }
2406  }
2407 
2408  if (coplanar) {
2409  /* excellent! - just make one big face */
2410  struct edgeuse *eu_tmp;
2411 
2412  /* put vertices in table */
2413  bu_ptbl_reset(&vert_tbl);
2414  for (i=0; i<loop_size; i++) {
2415  eu = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[i]);
2416  bu_ptbl_ins(&vert_tbl, (long *)&eu->vu_p->v_p);
2417  }
2418 
2419  /* make face */
2420  fu = nmg_cmface(s, (struct vertex ***)BU_PTBL_BASEADDR(&vert_tbl), loop_size);
2421 
2422  /* face geometry */
2423  if (nmg_loop_plane_area(BU_LIST_FIRST(loopuse, &fu->lu_hd), pl2) < 0.0) {
2424  bu_log("Failed planeeq\n");
2425  nmg_kfu(fu);
2426  } else {
2427  nmg_face_g(fu, pl2);
2428 
2429  /* Calculate face bounding box */
2430  nmg_face_bb(fu->f_p, tol);
2431  }
2432 
2433  /* now eliminate loop from table */
2434  eu_used = (struct edgeuse **)bu_calloc(loop_size, sizeof(struct edgeuse *), "edges used list");
2435  for (i=0; i<loop_size; i++)
2436  eu_used[i] = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[i]);
2437 
2438  for (i=0; i<loop_size; i++)
2439  bu_ptbl_rm(&eu_tbl, (long *)eu_used[i]);
2440 
2441  bu_free((char *)eu_used, "edge used list");
2442 
2443  /* may need to remove one more edgeuse from table */
2444  eu_tmp = nmg_find_e((*(struct vertex **)BU_PTBL_GET(&vert_tbl, 0)), (*(struct vertex **)BU_PTBL_GET(&vert_tbl, loop_size-1)), (struct shell *)NULL, (struct edge *)NULL);
2445  if (eu_tmp) {
2446  if (eu_tmp->radial_p != eu_tmp->eumate_p) {
2447  for (i=0; i<BU_PTBL_END(&eu_tbl); i++) {
2448  eu2 = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, i);
2449  if (EDGESADJ(eu2, eu_tmp)) {
2450  bu_ptbl_rm(&eu_tbl, (long *)eu2);
2451  break;
2452  }
2453  }
2454  }
2455  }
2456 
2457  /* set some flags to get us back to start of loop */
2458  found_face = 1;
2459  give_up_on_face = 1;
2460  }
2461 
2462  /* OK, so we have to do this one little-by-little */
2463  start_idx = -1;
2464  end_idx = -1;
2465 
2466  if (!found_face) {
2467 
2468  /* refresh the vertex list */
2469  (void)bu_ptbl_reset(&vert_tbl);
2470 
2471  start_idx++;
2472  end_idx = start_idx + 1;
2473  if (end_idx == loop_size)
2474  end_idx = 0;
2475 
2476  /* Get two edgeuses from the loop */
2477  eu1 = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[start_idx]);
2478  bu_ptbl_ins(&vert_tbl, (long *)&eu1->vu_p->v_p);
2479 
2480  eu2 = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[end_idx]);
2481  bu_ptbl_ins(&vert_tbl, (long *)&eu2->vu_p->v_p);
2482 
2483  edges_used = 2;
2484  /* if the edges are collinear, we can't make a face */
2485  while (bn_3pts_collinear(
2486  eu1->vu_p->v_p->vg_p->coord,
2487  eu2->vu_p->v_p->vg_p->coord,
2488  eu2->eumate_p->vu_p->v_p->vg_p->coord,
2489  tol) && edges_used < loop_size)
2490  {
2491  /* So, add another edge */
2492  end_idx++;
2493  if (end_idx == loop_size)
2494  end_idx = 0;
2495  eu1 = eu2;
2496  eu2 = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[end_idx]);
2497  bu_ptbl_ins(&vert_tbl, (long *)&eu2->vu_p->v_p);
2498  edges_used++;
2499  }
2500  }
2501 
2502  if (give_up_on_face) {
2503  loop_size = 0;
2504  start_loop = -1;
2505  break;
2506  }
2507 
2508  /* add last vertex to table */
2509  bu_ptbl_ins(&vert_tbl, (long *)&eu2->eumate_p->vu_p->v_p);
2510 
2511  /* save list of used edges to be removed later */
2512  eu_used = (struct edgeuse **)bu_calloc(edges_used, sizeof(struct edgeuse *), "edges used list");
2513  for (i=0; i<edges_used; i++)
2514  eu_used[i] = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[i]);
2515 
2516  /* make a face */
2517  fu = nmg_cmface(s, (struct vertex ***)BU_PTBL_BASEADDR(&vert_tbl), edges_used+1);
2518  /* out with the old, in with the new */
2519  for (i=0; i<edges_used; i++)
2520  bu_ptbl_rm(&eu_tbl, (long *)eu_used[i]);
2521 
2522  if (nmg_loop_plane_area(BU_LIST_FIRST(loopuse, &fu->lu_hd), pl2) < 0.0) {
2523  bu_log("Failed planeeq\n");
2524  nmg_kfu(fu);
2525  } else {
2526  nmg_face_g(fu, pl2);
2527 
2528  /* find the newly created edgeuse */
2529  found = 0;
2530  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
2531  NMG_CK_LOOPUSE(lu);
2532  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) ==
2534  continue;
2535  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
2536  NMG_CK_EDGEUSE(eu);
2537  if (eu->vu_p->v_p == (*(struct vertex **)BU_PTBL_GET(&vert_tbl, 0))
2538  && eu->eumate_p->vu_p->v_p == (*(struct vertex **)BU_PTBL_GET(&vert_tbl, edges_used)))
2539  {
2540  eu_new = eu;
2541  found = 1;
2542  break;
2543  } else if (eu->vu_p->v_p == (*(struct vertex **)BU_PTBL_GET(&vert_tbl, edges_used))
2544  && eu->eumate_p->vu_p->v_p == (*(struct vertex **)BU_PTBL_GET(&vert_tbl, 0)))
2545  {
2546  eu_new = eu->eumate_p;
2547  found = 1;
2548  break;
2549  }
2550 
2551  }
2552  if (found)
2553  break;
2554  }
2555 
2556  if (eu_new->radial_p == eu_new->eumate_p)
2557  bu_ptbl_ins(&eu_tbl, (long *)eu_new);
2558  else {
2559  struct edgeuse *eu_tmp;
2560 
2561  /* find third eu to be removed from eu_tbl */
2562  for (i=0; i<BU_PTBL_END(&eu_tbl); i++) {
2563  eu_tmp = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, i);
2564  if (eu_tmp->vu_p->v_p == eu_new->vu_p->v_p &&
2565  eu_tmp->eumate_p->vu_p->v_p == eu_new->eumate_p->vu_p->v_p)
2566  {
2567  bu_ptbl_rm(&eu_tbl, (long *)eu_tmp);
2568  break;
2569  }
2570  }
2571  }
2572  }
2573 
2574  bu_free((char *)eu_used, "edge used list");
2575 
2576  /* re-order loop */
2577  loop_size = 0;
2578  start_loop = -1;
2579  while (loop_size < 3 && start_loop < BU_PTBL_END(&eu_tbl)-3)
2580  order_tbl(&eu_tbl, ++start_loop, &idx, BU_PTBL_END(&eu_tbl), &loop_size);
2581  }
2582 
2583  if (give_up_on_face)
2584  continue;
2585 
2586  if (loop_size == 2) {
2587  bu_ptbl_reset(&vert_tbl);
2588  eu = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[0]);
2589  (void)bu_ptbl_ins(&vert_tbl, (long *)&eu->vu_p->v_p);
2590  eu = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[1]);
2591  (void)bu_ptbl_ins(&vert_tbl, (long *)&eu->vu_p->v_p);
2592  (void)bu_ptbl_ins(&vert_tbl, (long *)&eu->eumate_p->vu_p->v_p);
2593  if (!bn_3pts_collinear(
2594  ((*(struct vertex **)BU_PTBL_GET(&vert_tbl, 0)))->vg_p->coord,
2595  ((*(struct vertex **)BU_PTBL_GET(&vert_tbl, 1)))->vg_p->coord,
2596  ((*(struct vertex **)BU_PTBL_GET(&vert_tbl, 2)))->vg_p->coord,
2597  tol))
2598  {
2599  fu = nmg_cmface(s, (struct vertex ***)BU_PTBL_BASEADDR(&vert_tbl), 3);
2600  if (nmg_calc_face_g(fu)) {
2601  bu_log("Failed planeeq\n");
2602  nmg_kfu(fu);
2603  }
2604  } else
2605  bu_log("Not making face, edges are collinear!\n");
2606 
2607  loop_size = 0;
2608  continue;
2609  } else if (loop_size < 2) {
2610  loop_size = 0;
2611  continue;
2612  }
2613 
2614  /* if the last 3 vertices are collinear, then don't make the last face */
2615  bu_ptbl_reset(&vert_tbl);
2616  for (i=0; i<3; i++) {
2617  eu = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[i]);
2618  (void)bu_ptbl_ins(&vert_tbl, (long *)&eu->vu_p->v_p);
2619  }
2620 
2621  if (!bn_3pts_collinear(
2622  ((*(struct vertex **)BU_PTBL_GET(&vert_tbl, 0)))->vg_p->coord,
2623  ((*(struct vertex **)BU_PTBL_GET(&vert_tbl, 1)))->vg_p->coord,
2624  ((*(struct vertex **)BU_PTBL_GET(&vert_tbl, 2)))->vg_p->coord,
2625  tol))
2626  {
2627 
2628  /* Create last face from remaining 3 edges */
2629  fu = nmg_cmface(s, (struct vertex ***)BU_PTBL_BASEADDR(&vert_tbl), 3);
2630  if (nmg_calc_face_g(fu))
2631  bu_log("Failed planeeq\n");
2632 
2633  } else
2634  bu_log("Not making face, edges are collinear!\n");
2635 
2636  /* remove the last three edges from the table */
2637  {
2638  struct edgeuse *eu1, *eu2, *eu3;
2639 
2640  eu1 = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[0]);
2641  eu2 = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[1]);
2642  eu3 = (struct edgeuse *)BU_PTBL_GET(&eu_tbl, idx[2]);
2643  bu_ptbl_rm(&eu_tbl, (long *)eu1);
2644  bu_ptbl_rm(&eu_tbl, (long *)eu2);
2645  bu_ptbl_rm(&eu_tbl, (long *)eu3);
2646  }
2647  }
2648 
2649  /* Free up all the memory */
2650  bu_free((char *)idx, "idx");
2651  bu_ptbl_free(&eu_tbl);
2652  bu_ptbl_free(&vert_tbl);
2653 
2654 }
2655 
2656 
2657 /**
2658  * Duplicate a shell and return the new copy. New shell is in the same
2659  * region.
2660  *
2661  * The vertex geometry is copied from the source faces into
2662  * topologically distinct (new) vertex and vertex_g structs. They
2663  * will start out being geometrically coincident, but it is anticipated
2664  * that the caller will modify the geometry, e.g. as in an extrude
2665  * operation.
2666  *
2667  * NOTE: This routine creates a translation table that gives the
2668  * correspondence between old and new structures, the caller is
2669  * responsible for freeing this memory. Warning - NOT EVERY structure
2670  * is assigned a correspondence.
2671  */
2672 struct shell *
2673 nmg_dup_shell(struct shell *s, long int ***trans_tbl, const struct bn_tol *tol)
2674 {
2675 
2676  struct model *m;
2677  struct shell *new_s;
2678  struct faceuse *fu;
2679  struct loopuse *lu, *new_lu;
2680  struct edgeuse *eu;
2681  struct faceuse *new_fu;
2682  struct bu_ptbl faces;
2683  int tbl_size;
2684 
2685  if (RTG.NMG_debug & DEBUG_BASIC)
2686  bu_log("nmg_dup_shell(s = %p, trans_tbl = %p)\n", (void *)s, (void *)trans_tbl);
2687 
2688  NMG_CK_SHELL(s);
2689  BN_CK_TOL(tol);
2690 
2691  m = nmg_find_model((uint32_t *)s);
2692 
2693  /* create translation table double size to accommodate both copies */
2694  tbl_size = m->maxindex * 3;
2695  (*trans_tbl) = (long **)bu_calloc(tbl_size, sizeof(long *),
2696  "nmg_dup_shell trans_tbl");
2697 
2698  bu_ptbl_init(&faces, 64, " &faces ");
2699 
2700  new_s = nmg_msv(s->r_p);
2701  if (s->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2702  NMG_INDEX_ASSIGN((*trans_tbl), s, (long *)new_s);
2703  if (new_s->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2704  NMG_INDEX_ASSIGN((*trans_tbl), new_s, (long *)s);
2705 
2706  /* copy face uses */
2707  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
2708  NMG_CK_FACEUSE(fu);
2709  if (fu->orientation == OT_SAME) {
2710  new_fu = (struct faceuse *)NULL;
2711  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
2712  NMG_CK_LOOPUSE(lu);
2713  if (new_fu) {
2714  new_lu = nmg_dup_loop(lu, &new_fu->l.magic, (*trans_tbl));
2715  if (lu->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2716  NMG_INDEX_ASSIGN((*trans_tbl), lu, (long *)new_lu);
2717  if (new_lu->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2718  NMG_INDEX_ASSIGN((*trans_tbl), new_lu, (long *)lu);
2719  } else {
2720  new_lu = nmg_dup_loop(lu, &new_s->l.magic, (*trans_tbl));
2721  if (new_lu->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2722  NMG_INDEX_ASSIGN((*trans_tbl), lu, (long *)new_lu);
2723  if (lu->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2724  NMG_INDEX_ASSIGN((*trans_tbl), new_lu, (long *)lu);
2725  new_fu = nmg_mf(new_lu);
2726  if (lu->orientation == OT_OPPOSITE) {
2727  /* nmg_mf forces loops to OT_SAME */
2728  new_lu->orientation = OT_OPPOSITE;
2729  new_lu->lumate_p->orientation = OT_OPPOSITE;
2730  }
2731  if (fu->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2732  NMG_INDEX_ASSIGN((*trans_tbl), fu, (long *)new_fu);
2733  if (new_fu->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2734  NMG_INDEX_ASSIGN((*trans_tbl), new_fu, (long *)fu);
2735  if (fu->fumate_p->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2736  NMG_INDEX_ASSIGN((*trans_tbl), fu->fumate_p, (long *)new_fu->fumate_p);
2737  if (new_fu->fumate_p->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2738  NMG_INDEX_ASSIGN((*trans_tbl), new_fu->fumate_p, (long *)fu->fumate_p);
2739  if (fu->f_p->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2740  NMG_INDEX_ASSIGN((*trans_tbl), fu->f_p, (long *)new_fu->f_p);
2741  if (new_fu->f_p->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2742  NMG_INDEX_ASSIGN((*trans_tbl), new_fu->f_p, (long *)fu->f_p);
2743  }
2744  }
2745 
2746  /* make sure it's not a face without a loop before proceeding */
2747  if (!new_fu)
2748  continue;
2749 
2750  if (fu->f_p->g.plane_p) {
2751  /* Do it this way if you expect to change the normals */
2752  plane_t n;
2753  NMG_GET_FU_PLANE(n, fu);
2754  nmg_face_g(new_fu, n);
2755 
2756  /* XXX Perhaps this should be new_fu->f_p->g.plane_p ? */
2757  if (fu->f_p->g.plane_p->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2758  NMG_INDEX_ASSIGN((*trans_tbl), fu->f_p->g.plane_p, (long *)new_fu->f_p->g.plane_p);
2759  if (new_fu->f_p->g.plane_p->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2760  NMG_INDEX_ASSIGN((*trans_tbl), new_fu->f_p->g.plane_p, (long *)fu->f_p->g.plane_p);
2761  }
2762  new_fu->orientation = fu->orientation;
2763  new_fu->fumate_p->orientation = fu->fumate_p->orientation;
2764  bu_ptbl_ins(&faces, (long *)new_fu);
2765  }
2766  }
2767 
2768  /* glue new faces */
2769  nmg_gluefaces((struct faceuse **)BU_PTBL_BASEADDR(&faces), BU_PTBL_END(&faces), tol);
2770  bu_ptbl_free(&faces);
2771 
2772  /* copy wire loops */
2773  for (BU_LIST_FOR (lu, loopuse, &s->lu_hd)) {
2774  NMG_CK_LOOPUSE(lu);
2775  new_lu = nmg_dup_loop(lu, &new_s->l.magic, (*trans_tbl));
2776  if (lu->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2777  NMG_INDEX_ASSIGN((*trans_tbl), lu, (long *)new_lu);
2778  if (new_lu->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2779  NMG_INDEX_ASSIGN((*trans_tbl), new_lu, (long *)lu);
2780  }
2781 
2782  /* copy wire edges */
2783  for (BU_LIST_FOR (eu, edgeuse, &s->eu_hd)) {
2784  struct vertex *old_v1, *old_v2, *new_v1, *new_v2;
2785  struct edgeuse *new_eu;
2786 
2787  NMG_CK_EDGEUSE(eu);
2788  NMG_CK_VERTEXUSE(eu->vu_p);
2789  NMG_CK_VERTEX(eu->vu_p->v_p);
2790  NMG_CK_EDGEUSE(eu->eumate_p);
2791  NMG_CK_VERTEXUSE(eu->eumate_p->vu_p);
2792  NMG_CK_VERTEX(eu->eumate_p->vu_p->v_p);
2793 
2794  old_v1 = eu->vu_p->v_p;
2795  new_v1 = NMG_INDEX_GETP(vertex, (*trans_tbl), old_v1);
2796  old_v2 = eu->eumate_p->vu_p->v_p;
2797  new_v2 = NMG_INDEX_GETP(vertex, (*trans_tbl), old_v2);
2798 
2799  /* make the wire edge */
2800  new_eu = nmg_me(new_v1, new_v2, new_s);
2801  if (eu->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2802  NMG_INDEX_ASSIGN((*trans_tbl), eu, (long *)new_eu);
2803  if (new_eu->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2804  NMG_INDEX_ASSIGN((*trans_tbl), new_eu, (long *)eu);
2805 
2806  new_v1 = new_eu->vu_p->v_p;
2807  if (old_v1->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2808  NMG_INDEX_ASSIGN((*trans_tbl), old_v1, (long *)new_v1);
2809  if (new_v1->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2810  NMG_INDEX_ASSIGN((*trans_tbl), new_v1, (long *)old_v1);
2811  if (!new_v1->vg_p) {
2812  nmg_vertex_gv(new_v1, old_v1->vg_p->coord);
2813  if (old_v1->vg_p->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2814  NMG_INDEX_ASSIGN((*trans_tbl), old_v1->vg_p, (long *)new_v1->vg_p);
2815  if (new_v1->vg_p->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2816  NMG_INDEX_ASSIGN((*trans_tbl), new_v1->vg_p, (long *)old_v1->vg_p);
2817  }
2818 
2819  new_v2 = new_eu->eumate_p->vu_p->v_p;
2820  if (old_v2->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2821  NMG_INDEX_ASSIGN((*trans_tbl), old_v2, (long *)new_v2);
2822  if (new_v2->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2823  NMG_INDEX_ASSIGN((*trans_tbl), new_v2, (long *)old_v2);
2824  if (!new_v2->vg_p) {
2825  nmg_vertex_gv(new_v2, old_v2->vg_p->coord);
2826  if (old_v2->vg_p->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2827  NMG_INDEX_ASSIGN((*trans_tbl), old_v2->vg_p, (long *)new_v2->vg_p);
2828  if (new_v2->vg_p->index >= tbl_size) bu_bomb("nmg_dup_shell: trans table exceeded\n");
2829  NMG_INDEX_ASSIGN((*trans_tbl), new_v2->vg_p, (long *)old_v2->vg_p);
2830  }
2831 
2832  }
2833 
2834 
2835  return new_s;
2836 }
2837 
2838 
2839 /* Routines to use the bu_ptbl structures as a stack of edgeuse structures*/
2840 #define NMG_PUSH(_ptr, _stack) bu_ptbl_ins(_stack, (long *) _ptr);
2841 
2842 struct edgeuse
2843 *nmg_pop_eu(struct bu_ptbl *stack)
2844 {
2845  struct edgeuse *eu;
2846 
2847  /* return a NULL if stack is empty */
2848  if (BU_PTBL_END(stack) == 0)
2849  return (struct edgeuse *)NULL;
2850 
2851  /* get last edgeuse on the stack */
2852  eu = (struct edgeuse *)BU_PTBL_GET(stack, BU_PTBL_END(stack)-1);
2853 
2854  /* remove that edgeuse from the stack */
2855  bu_ptbl_rm(stack, (long *)eu);
2856 
2857  return eu;
2858 }
2859 
2860 
2861 void
2862 nmg_reverse_radials(struct faceuse *fu, const struct bn_tol *tol)
2863 {
2864  struct loopuse *lu;
2865 
2866  if (RTG.NMG_debug & DEBUG_BASIC)
2867  bu_log("nmg_reverse_radials(fu = %p)\n", (void *)fu);
2868 
2869  NMG_CK_FACEUSE(fu);
2870  BN_CK_TOL(tol);
2871 
2872  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
2873  struct edgeuse *eu;
2874  struct edgeuse *eu_radial;
2875  struct edgeuse *eumate;
2876  struct edgeuse *eumate_radial;
2877 
2878  NMG_CK_LOOPUSE(lu);
2879 
2880  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
2881  continue;
2882 
2883  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
2884  eu_radial = eu->radial_p;
2885  eumate = eu->eumate_p;
2886  eumate_radial = eumate->radial_p;
2887 
2888  /* if no radials continue to next edgeuse in loop */
2889  if (eu_radial == eumate)
2890  continue;
2891 
2892  /* Note: there is no problem if this loop is radial to another
2893  * loop in the same face, the radials will get switched as we process
2894  * the first loop, then switched back as we process the second
2895  */
2896 
2897  eu_radial->radial_p = eumate;
2898  eu->radial_p = eumate_radial;
2899  eumate_radial->radial_p = eu;
2900  eumate->radial_p = eu_radial;
2901  }
2902  }
2903 }
2904 
2905 
2906 /**
2907  * This routine calls "nmg_reverse_face" and also makes the radial
2908  * pointers connect faces of like orientation (i.e., OT_SAME to
2909  * OT_SAME and OT_OPPOSITE to OT_OPPOSITE).
2910  *
2911  * XXX Don't use this, use nmg_s_radial_harmonize() at the right time.
2912  */
2913 void
2914 nmg_reverse_face_and_radials(struct faceuse *fu, const struct bn_tol *tol)
2915 {
2916  struct loopuse *lu;
2917 
2918  if (RTG.NMG_debug & DEBUG_BASIC)
2919  bu_log("nmg_reverse_face_and_radials(fu = %p)\n", (void *)fu);
2920 
2921  NMG_CK_FACEUSE(fu);
2922  BN_CK_TOL(tol);
2923 
2924  /* reverse face */
2925  nmg_reverse_face(fu);
2926 
2927  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
2928  struct edgeuse *eu;
2929  struct edgeuse *eu_radial;
2930  struct edgeuse *eumate;
2931  struct edgeuse *eumate_radial;
2932 
2933  NMG_CK_LOOPUSE(lu);
2934 
2935  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
2936  continue;
2937 
2938  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
2939  eu_radial = eu->radial_p;
2940  eumate = eu->eumate_p;
2941  eumate_radial = eumate->radial_p;
2942 
2943  /* if no radials continue to next edgeuse in loop */
2944  if (eu_radial == eumate)
2945  continue;
2946 
2947  /* Note: there is no problem if this loop is radial to another
2948  * loop in the same face, the radials will get switched as we process
2949  * the first loop, then switched back as we process the second
2950  */
2951 
2952  eu_radial->radial_p = eumate;
2953  eu->radial_p = eumate_radial;
2954  eumate_radial->radial_p = eu;
2955  eumate->radial_p = eu_radial;
2956  }
2957  }
2958 }
2959 
2960 
2961 /**
2962  * determines if the shell is a void shell or an exterior shell by
2963  * finding the topmost face (in some direction) and looking at the
2964  * component of the OT_SAME faceuse normal in that direction.
2965  *
2966  * returns:
2967  * 0 - shell is exterior shell
2968  * 1 - shell is a void shell
2969  * -1 - cannot determine
2970  *
2971  * It is expected that this shell is the result of
2972  * nmg_decompose_shells().
2973  */
2974 int
2975 nmg_shell_is_void(const struct shell *s)
2976 {
2977  struct model *m;
2978  struct face *f;
2979  struct faceuse *fu;
2980  vect_t normal;
2981  int dir;
2982  long *flags;
2983 
2984  NMG_CK_SHELL(s);
2985 
2986  m = nmg_find_model(&s->l.magic);
2987  NMG_CK_MODEL(m);
2988 
2989  flags = (long *)bu_calloc(m->maxindex, sizeof(long), "nmg_shell_is_void: flags ");
2990 
2991  f = nmg_find_top_face(s, &dir, flags);
2992 
2993  bu_free((char *)flags, "nmg_shell_is_void: flags");
2994 
2995  if (f == (struct face *)NULL)
2996  return -1;
2997 
2998  NMG_CK_FACE(f);
2999  NMG_CK_FACE_G_PLANE(f->g.plane_p);
3000  fu = f->fu_p;
3001  NMG_CK_FACEUSE(fu);
3002 
3003  if (fu->orientation != OT_SAME)
3004  fu = fu->fumate_p;
3005  if (fu->orientation != OT_SAME)
3006  return -1;
3007 
3008  NMG_GET_FU_NORMAL(normal, fu);
3009 
3010  if (ZERO(normal[dir]))
3011  return -1;
3012  if (normal[dir] < 0.0)
3013  return 1;
3014  else
3015  return 0;
3016 }
3017 
3018 
3019 /**
3020  * This routine expects "fu_in" to have a correctly oriented normal.
3021  * It then checks all faceuses in the same shell it can reach via
3022  * radial structures, and reverses faces and modifies radial
3023  * structures as needed to result in a consistent NMG shell
3024  * structure. The "flags" variable is a translation table for the
3025  * model, and as each face is checked, its flag is set. Faces with
3026  * flags that have already been set will not be checked by this
3027  * routine.
3028  */
3029 void
3030 nmg_propagate_normals(struct faceuse *fu_in, long int *flags, const struct bn_tol *tol)
3031 {
3032  struct bu_ptbl stack;
3033  struct loopuse *lu;
3034  struct edgeuse *eu, *eu1;
3035  struct faceuse *fu;
3036 
3037  if (RTG.NMG_debug & DEBUG_BASIC)
3038  bu_log("nmg_propagate_normals(fu_in = %p, flags = %p)\n", (void *)fu_in, (void *)flags);
3039 
3040  NMG_CK_FACEUSE(fu_in);
3041  BN_CK_TOL(tol);
3042 
3043  fu = fu_in;
3044  if (fu->orientation != OT_SAME)
3045  fu = fu->fumate_p;
3046  if (fu->orientation != OT_SAME) {
3047  bu_log("nmg_propagate_normals: Could not find OT_SAME orientation of faceuse %p\n", (void *)fu_in);
3048  return;
3049  }
3050 
3051  /* set flag for this face since we know this one is OK */
3052  NMG_INDEX_SET(flags, fu->f_p);
3053 
3054  /* Use the ptbl structure as a stack */
3055  bu_ptbl_init(&stack, 64, " &stack ");
3056 
3057  /* push all edgeuses of "fu" onto the stack */
3058  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
3059  NMG_CK_LOOPUSE(lu);
3060  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
3061  continue;
3062  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
3063  /* don't push free edges on the stack */
3064  if (eu->radial_p->eumate_p != eu)
3065  NMG_PUSH(eu, &stack);
3066  }
3067  }
3068 
3069  /* now pop edgeuses from stack, go to radial face, fix its normal,
3070  * and push its edgeuses onto the stack */
3071 
3072  while ((eu1 = nmg_pop_eu(&stack)) != (struct edgeuse *)NULL) {
3073  /* eu1 is an edgeuse on an OT_SAME face, so its radial
3074  * should be in an OT_SAME also */
3075 
3076  NMG_CK_EDGEUSE(eu1);
3077 
3078  eu = eu1->radial_p;
3079  NMG_CK_EDGEUSE(eu);
3080 
3081  /* find the face that contains this edgeuse */
3082  fu = nmg_find_fu_of_eu(eu);
3083  if (!fu)
3084  continue;
3085 
3086  NMG_CK_FACEUSE(fu);
3087 
3088  /* if this face has already been processed, skip it */
3089  if (NMG_INDEX_TEST_AND_SET(flags, fu->f_p)) {
3090  if (RTG.NMG_debug & DEBUG_BASIC)
3091  bu_log("nmg_propagate_normals: checking fu %p, eu = %p, eu1 = %p\n",
3092  (void *)fu, (void *)eu, (void *)eu1);
3093 
3094  if (fu->orientation == OT_SAME) {
3095  if (eu1->vu_p->v_p == eu->vu_p->v_p &&
3096  eu1->eumate_p->vu_p->v_p == eu->eumate_p->vu_p->v_p) {
3097  /* edge direction is wrong, need to reverse
3098  * both the face and the radials
3099  */
3101  }
3102  /* else - this is the way it should be */
3103  } else if (fu->orientation == OT_OPPOSITE) {
3104  if (eu1->vu_p->v_p == eu->vu_p->v_p &&
3105  eu1->eumate_p->vu_p->v_p == eu->eumate_p->vu_p->v_p)
3106  {
3107  /* just swap radial pointer around */
3108  nmg_reverse_radials(fu, tol);
3109  } else {
3110  /* just reverse the face */
3111  nmg_reverse_face(fu);
3112  }
3113  } else {
3114  bu_log("nmg_propagate_normals: found an unoriented face!\n");
3115  nmg_pr_fu_briefly(fu, "");
3116  bu_bomb("nmg_propagate_normals: found an unoriented face!\n");
3117  }
3118 
3119  /* make sure we are dealing with an OT_SAME faceuse */
3120  if (fu->orientation != OT_SAME)
3121  fu = fu->fumate_p;
3122 
3123  /* push all edgeuses of "fu" onto the stack */
3124  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
3125  NMG_CK_LOOPUSE(lu);
3126  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
3127  continue;
3128  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
3129  /* don't push free edges on the stack */
3130  if (eu->radial_p->eumate_p != eu)
3131  NMG_PUSH(eu, &stack);
3132  }
3133  }
3134  }
3135  }
3136 
3137  /* free the stack */
3138  bu_ptbl_free(&stack);
3139 }
3140 
3141 
3142 /**
3143  * looks for edges that have uses in more than one shell in region.
3144  * creates new edges so that there is no sharing of edges among shells
3145  */
3146 HIDDEN void
3147 nmg_disconnect_shells(struct nmgregion *r)
3148 {
3149  struct shell *s1;
3150  struct bu_ptbl edges;
3151  int count=0;
3152  int i;
3153 
3154  NMG_CK_REGION(r);
3155 
3156  /* count number of shells in this region */
3157  for (BU_LIST_FOR (s1, shell, &r->s_hd))
3158  count++;
3159 
3160  /* if there is less than two shells, nothing to do */
3161  if (count < 2)
3162  return;
3163 
3164  /* get a list of all edges in this region */
3165  nmg_edge_tabulate(&edges, &r->l.magic);
3166 
3167  /* look at every edge in region */
3168  for (i=0; i<BU_PTBL_END(&edges); i++) {
3169  struct edge *e;
3170  struct edgeuse *eu;
3171  struct edgeuse *eu1;
3172  struct shell *needs_disconnect=(struct shell *)NULL;
3173 
3174  e = (struct edge *)BU_PTBL_GET(&edges, i);
3175  NMG_CK_EDGE(e);
3176 
3177  /* find shell of a use of this edge */
3178  eu = e->eu_p;
3179  s1 = nmg_find_s_of_eu(eu);
3180 
3181  /* check if any other use of this edge is from a different shell */
3182  eu1 = eu->eumate_p->radial_p;
3183  while (eu1 != eu && eu1 != eu->eumate_p) {
3184  struct shell *s2;
3185 
3186  if ((s2 = nmg_find_s_of_eu(eu1)) != s1) {
3187  needs_disconnect = s2;
3188  break;
3189  }
3190  eu1 = eu1->eumate_p->radial_p;
3191  }
3192 
3193  while (needs_disconnect) {
3194  struct edgeuse *eu2;
3195  struct edgeuse *start_eu;
3196  struct edgeuse *last;
3197  int last_orientation;
3198 
3199  /* disconnect first use of this edge in shell 'needs_disconnect' */
3200  start_eu = eu1->radial_p;
3201  nmg_unglueedge(eu1);
3202  last = eu1;
3203  last_orientation = (nmg_find_fu_of_eu(eu1))->orientation;
3204 
3205  /* now disconnect all other uses, reconnecting them to eu1 as we go */
3206  while (nmg_find_s_of_eu(start_eu) == needs_disconnect)
3207  start_eu = start_eu->eumate_p->radial_p;
3208  eu2 = start_eu;
3209  do {
3210  struct edgeuse *next_eu;
3211  struct faceuse *fu2;
3212 
3213  /* find uses in 'needs_disconnect' shell */
3214  next_eu = eu2->eumate_p->radial_p;
3215  if (nmg_find_s_of_eu(eu2) == needs_disconnect) {
3216 
3217  /* disconnect this use */
3218  nmg_unglueedge(eu2);
3219  fu2 = nmg_find_fu_of_eu(eu2);
3220 
3221  /* reconnect it to 'needs_disconnect' shell */
3222  if (fu2->orientation == last_orientation) {
3223  nmg_je(last, eu2);
3224  last = eu2->eumate_p;
3225  last_orientation = fu2->fumate_p->orientation;
3226  } else {
3227  nmg_je(last, eu2->eumate_p);
3228  last_orientation = fu2->orientation;
3229  last = eu2;
3230  }
3231  }
3232  eu2 = next_eu;
3233  } while (eu2 != start_eu && eu2->eumate_p != start_eu);
3234 
3235  /* now check remaining uses */
3236  eu1 = eu->eumate_p->radial_p;
3237  needs_disconnect = (struct shell *)NULL;
3238  while (eu1 != eu && eu1 != eu->eumate_p) {
3239  struct shell *s2;
3240 
3241  if ((s2 = nmg_find_s_of_eu(eu1)) != s1) {
3242  needs_disconnect = s2;
3243  break;
3244  }
3245  eu1 = eu1->eumate_p->radial_p;
3246  }
3247  }
3248 
3249  }
3250 }
3251 
3252 
3253 /**
3254  * looks for radially connected faceuses that have disagreeing
3255  * orientations. if such a condition is found, the radial pointers
3256  * are rearranged to make the radial faces agree in orientation.
3257  */
3258 void
3260 {
3261  struct faceuse *fu;
3262 
3263  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
3264  struct loopuse *lu;
3265 
3266  if (fu->orientation != OT_SAME)
3267  continue;
3268 
3269  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
3270  struct edgeuse *eu;
3271 
3272  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
3273  continue;
3274 
3275  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
3276  struct faceuse *fu2;
3277  struct faceuse *fu3;
3278  struct edgeuse *eu2;
3279 
3280  eu2 = eu->radial_p;
3281 
3282  if (eu2 == eu->eumate_p)
3283  continue;
3284 
3285  fu2 = nmg_find_fu_of_eu(eu2);
3286  if (fu2->orientation == OT_SAME)
3287  continue;
3288 
3289  fu3 = nmg_find_fu_of_eu(eu->eumate_p->radial_p);
3290  if (fu3->orientation == OT_OPPOSITE)
3291  continue;
3292  }
3293  }
3294  }
3295 }
3296 
3297 
3298 /**
3299  * Routine to set all OT_SAME faceuse normals to outward direction.
3300  * Assumes that there are no other shells sharing edges with this one.
3301  */
3302 void
3303 nmg_fix_decomposed_shell_normals(struct shell *s, const struct bn_tol *tol)
3304 {
3305  struct model *m;
3306  struct face *f_top;
3307  struct faceuse *fu;
3308  vect_t normal;
3309  int dir;
3310  int missed_faces;
3311  long *flags;
3312 
3313  NMG_CK_SHELL(s);
3314  BN_CK_TOL(tol);
3315 
3316  m = s->r_p->m_p;
3317  flags = (long *)bu_calloc(m->maxindex, sizeof(long), "nmg_fix_decomposed_shell_normals: flags");
3318 
3319 missed:
3320  /* find the top face */
3321  f_top = nmg_find_top_face(s, &dir, flags);
3322  if (f_top == (struct face *)NULL) {
3323  bu_log("nmg_fix_decomposed_shell_normals: Could not get a top face from nmg_find_top_face()\n");
3324  bu_log("\tWARNING: continuing without fixing normals!\n");
3325  bu_free((char *)flags, "nmg_fix_decomposed_shell_normals: flags");
3326  return;
3327  }
3328  if (*f_top->g.magic_p != NMG_FACE_G_PLANE_MAGIC) {
3329  NMG_INDEX_SET(flags, f_top);
3330  goto missed;
3331  }
3332 
3333  if (NMG_INDEX_TEST(flags, f_top))
3334  bu_log(" nmg_find_top_face returned a flagged face %p\n", (void *)f_top);
3335 
3336  NMG_CK_FACE(f_top);
3337  fu = f_top->fu_p;
3338  NMG_CK_FACEUSE(fu);
3339  if (fu->orientation != OT_SAME)
3340  fu = fu->fumate_p;
3341  if (fu->orientation != OT_SAME) {
3342  bu_log("nmg_fix_decomposed_shell_normals: no OT_SAME use of top face\n");
3343  bu_free((char *)flags, "nmg_fix_decomposed_shell_normals: flags");
3344  return;
3345  }
3346  NMG_GET_FU_NORMAL(normal, fu);
3347 
3348  if (RTG.NMG_debug & DEBUG_BASIC) {
3349  bu_log("\tnmg_fix_decomposed_shell_normals: top face is %p in %d direction, OT_SAME use is %p\n",
3350  (void *)f_top, dir, (void *)fu);
3351  bu_log("\toutward normal = (%g %g %g)\n", V3ARGS(normal));
3352  }
3353 
3354  /* f_top is the topmost face (in the "dir" direction), so its OT_SAME use should have a
3355  * normal with component in the "dir" direction */
3356  if (normal[dir] < 0.0) {
3357  if (RTG.NMG_debug & DEBUG_BASIC)
3358  bu_log("nmg_fix_decomposed_shell_normals: reversing fu %p\n", (void *)fu);
3359 
3361  }
3362 
3363  /* get OT_SAME use of top face */
3364  fu = f_top->fu_p;
3365  if (fu->orientation != OT_SAME)
3366  fu = fu->fumate_p;
3367 
3368  NMG_CK_FACEUSE(fu);
3369 
3370  /* fu is now known to be a correctly oriented faceuse,
3371  * propagate this throughout the shell, face by face, by
3372  * traversing the shell using the radial edge structure */
3373 
3374  nmg_propagate_normals(fu, flags, tol);
3375 
3376  if (RTG.NMG_debug & DEBUG_BASIC) {
3377  vect_t new_norm;
3378 
3379  NMG_GET_FU_NORMAL(new_norm, fu);
3380  bu_log("nmg_fix_decomposed_shell_normals: After propagation top faceuse normal is (%g %g %g)\n",
3381  V3ARGS(new_norm));
3382  }
3383 
3384  /* check if all the faces have been processed */
3385  missed_faces = 0;
3386  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
3387  NMG_CK_FACEUSE(fu);
3388  if (fu->orientation == OT_SAME) {
3389  if (!NMG_INDEX_TEST(flags, fu->f_p))
3390  missed_faces++;
3391  }
3392  }
3393 
3394  if (missed_faces) {
3395  bu_log("nmg_fix_decomposed_shell_normals: missed %d faces in shell %p (was it decomposed?)\n",
3396  missed_faces, (void *)s);
3397  bu_bomb("nmg_fix_decomposed_shell_normals: missed faces in shell (was it decomposed?)\n");
3398  }
3399 
3400  bu_free((char *)flags, "flags");
3401 }
3402 
3403 
3404 /*
3405  * Creates a new model from an existing nmgregion. Will refuse to
3406  * create new model if the passed nmgregion has children with uses in
3407  * another nmgregion.
3408  *
3409  * The reindex flag indicates if the new model should be reindexed.
3410  * If non-zero, nmg_m_reindex() is called. If zero, the maxindex
3411  * field of the new model is filled by finding the largest index value
3412  * in the new model and adding one. This is useful if the indices in
3413  * the region are being used for table look-up.
3414  */
3415 struct model *
3416 nmg_mk_model_from_region(struct nmgregion *r, int reindex)
3417 {
3418  struct model *m;
3419  struct bu_ptbl tbl;
3420  int i;
3421  int other_uses=0;
3422 
3423  NMG_CK_REGION(r);
3424 
3425 
3426  /* check if anything in this region has uses in another region */
3427  nmg_vertex_tabulate(&tbl, &r->l.magic);
3428 
3429  for (i=0; i<BU_PTBL_END(&tbl); i++) {
3430  struct vertex *v;
3431  struct vertexuse *vu;
3432 
3433  v = (struct vertex *)BU_PTBL_GET(&tbl, i);
3434  NMG_CK_VERTEX(v);
3435 
3436  for (BU_LIST_FOR (vu, vertexuse, &v->vu_hd)) {
3437  if ((nmg_find_s_of_vu(vu))->r_p != r) {
3438  bu_log("vertexuse %p (v=%p) at (%g %g %g) has use in another region\n",
3439  (void *)vu, (void *)v, V3ARGS(v->vg_p->coord));
3440  other_uses++;
3441  }
3442  }
3443  }
3444  bu_ptbl_free(&tbl);
3445 
3446  nmg_edge_tabulate(&tbl, &r->l.magic);
3447  for (i=0; i<BU_PTBL_END(&tbl); i++) {
3448  struct edge *e;
3449  struct edgeuse *eu;
3450  struct edgeuse *eu1;
3451 
3452  e = (struct edge *)BU_PTBL_GET(&tbl, i);
3453  NMG_CK_EDGE(e);
3454 
3455  eu = e->eu_p;
3456  NMG_CK_EDGEUSE(eu);
3457 
3458  eu1 = eu->radial_p->eumate_p;
3459  while (eu1 != eu) {
3460  if ((nmg_find_s_of_eu(eu1))->r_p != r) {
3461  bu_log("edgeuse %p (e=%p) at (%g %g %g)<->(%g %g %g0 has use in another region\n",
3462  (void *)eu, (void *)e, V3ARGS(eu->vu_p->v_p->vg_p->coord),
3463  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord));
3464  other_uses++;
3465  }
3466  eu1 = eu1->radial_p->eumate_p;
3467  }
3468  }
3469  bu_ptbl_free(&tbl);
3470 
3471  if (other_uses) {
3472  return (struct model *)NULL;
3473  }
3474 
3475  m = nmg_mm();
3476 
3477  BU_LIST_DEQUEUE(&r->l);
3478 
3479  BU_LIST_INSERT(&m->r_hd, &r->l);
3480 
3481  r->m_p = m;
3482 
3483  if (reindex)
3484  nmg_m_reindex(m, 0);
3485  else
3486  m->maxindex = nmg_find_max_index(m) + 1;
3487 
3488  return m;
3489 }
3490 
3491 
3492 /**
3493  * Routine to set faceuse normals to correct direction.
3494  *
3495  * Method:
3496  * 1. Make a copy of the shell in another model.
3497  * 2. Decompose the copy into constituent shells.
3498  * 3. Set all normals of constituent shells outward (no void shells)
3499  * 4. Use nmg_classify_s_vs_s() for every pair of shells.
3500  * 5. Mark any shell that is inside an odd number of other shells as a void shell.
3501  * 6. Compare original faceuses with the results in the copy and adjust normals as needed.
3502  * 7. Destroy the copy model.
3503  */
3504 void
3505 nmg_fix_normals(struct shell *s_orig, const struct bn_tol *tol)
3506 {
3507  struct model *tmp_m;
3508  struct model *m;
3509  struct shell *dup_s;
3510  struct shell *s1;
3511  struct nmgregion *tmp_r;
3512  struct faceuse *fu;
3513  struct bu_ptbl reverse;
3514  int shell_count;
3515  long **trans_tbl;
3516 
3517  if (RTG.NMG_debug & DEBUG_BASIC)
3518  bu_log("nmg_fix_normals(s = %p)\n", (void *)s_orig);
3519 
3520  NMG_CK_SHELL(s_orig);
3521  BN_CK_TOL(tol);
3522 
3523  /* Currently we can only fix normals for planar faces
3524  * check that there are no TNURB faces
3525  */
3526  for (BU_LIST_FOR (fu, faceuse, &s_orig->fu_hd)) {
3527  struct face *f;
3528 
3529  NMG_CK_FACEUSE(fu);
3530 
3531  if (fu->orientation != OT_SAME)
3532  continue;
3533 
3534  f = fu->f_p;
3535 
3536  if (!f->g.magic_p) {
3537  bu_log("nmg_fix_normals failed, found a face with no geometry (%p)\n", (void *)f);
3538  return;
3539  }
3540 
3541  if (*f->g.magic_p != NMG_FACE_G_PLANE_MAGIC) {
3542  bu_log("nmg_fix_normals: non-planar face found (%p)\n", (void *)f);
3543  bu_log(" cannot fix normals\n");
3544  return;
3545  }
3546  }
3547 
3548  m = s_orig->r_p->m_p;
3549 
3550  /* make a temporary nmgregion for us to work in */
3551  tmp_r = nmg_mrsv(m);
3552 
3553  /* get rid of the automatically created shell */
3554  (void)nmg_ks(BU_LIST_FIRST(shell, &tmp_r->s_hd));
3555 
3556  /* make a copy of the shell of interest */
3557  dup_s = nmg_dup_shell(s_orig, &trans_tbl, tol);
3558 
3559  /* move the copy to our work area */
3560  nmg_mv_shell_to_region(dup_s, tmp_r);
3561 
3562  /* move duplicate shell to another model */
3563  tmp_m = nmg_mk_model_from_region(tmp_r, 0); /* don't reindex, We need the old indices */
3564  nmg_rebound(tmp_m, tol);
3565 
3566  /* decompose the shell */
3567  shell_count = nmg_decompose_shell(dup_s, tol);
3568 
3569  if (shell_count == 1) {
3570  /* just one shell, so fix it and return */
3571  (void)nmg_km(tmp_m);
3572  bu_free((char *)trans_tbl, "translate table");
3574  nmg_fix_decomposed_shell_normals(s_orig, tol);
3575  return;
3576  }
3577 
3578  /* make sure the shells don't share any edges */
3579  nmg_disconnect_shells(tmp_r);
3580 
3581  /* Make sure all OT_SAME faceuses are radial to OT_SAME faceuses */
3582  for (BU_LIST_FOR (s1, shell, &tmp_r->s_hd))
3584 
3585  /* Decomposed into more than one shell.
3586  * Need to check for inner void shells.
3587  * Start by making all the shells look like solids (no voids).
3588  */
3589  for (BU_LIST_FOR (s1, shell, &tmp_r->s_hd))
3591 
3592  /* initialize a list of shells to be reversed */
3593  bu_ptbl_init(&reverse, 8, "Ptbl for nmg_fix_normals");
3594 
3595  /* now check which shells are inside others */
3596  for (BU_LIST_FOR (s1, shell, &tmp_r->s_hd)) {
3597  struct shell *s2;
3598  int inner_count=0;
3599  int stop = 0;
3600 
3601  for (BU_LIST_FOR (s2, shell, &tmp_r->s_hd)) {
3602  int nmg_class;
3603 
3604  if (s1 == s2)
3605  continue;
3606 
3607  nmg_class = nmg_classify_s_vs_s(s1, s2, tol);
3608  if (nmg_class == NMG_CLASS_AinB)
3609  inner_count++;
3610  else if (nmg_class == NMG_CLASS_Unknown) {
3611  bu_log("nmg_fix_normals: nmg_classify_s_vs_s() failed for shells %p and %p\n",
3612  (void *)s1, (void *)s2);
3613  bu_log(" Continuing anyway (shell is likely to have incorrectly oriented normals)\n");
3614  stop = 1;
3615  break;
3616  }
3617  }
3618 
3619  if (inner_count % 2) {
3620  /* shell s1 is inside an odd number of shells, so it must be a void */
3621  bu_ptbl_ins(&reverse, (long *)s1);
3622  }
3623  if (stop) {
3624  break;
3625  }
3626  }
3627 
3628  /* now set faces in original shell to match our calculations */
3630 
3631  for (BU_LIST_FOR (s1, shell, &tmp_r->s_hd)) {
3632  int reversed;
3633 
3634  if (bu_ptbl_locate(&reverse, (long *)s1) == (-1))
3635  reversed = 0;
3636  else
3637  reversed = 1;
3638 
3639  for (BU_LIST_FOR (fu, faceuse, &s1->fu_hd)) {
3640  struct faceuse *fu_in_s;
3641  vect_t normal;
3642  vect_t normal_in_s;
3643 
3644  if (fu->orientation != OT_SAME)
3645  continue;
3646 
3647  fu_in_s = NMG_INDEX_GETP(faceuse, trans_tbl, fu);
3648  if (!fu_in_s) {
3649  bu_log("fu %p does not have correspondence in original shell\n", (void *)fu);
3650  nmg_pr_fu_briefly(fu, "");
3651  continue;
3652  }
3653  if (fu_in_s->orientation != OT_SAME)
3654  fu_in_s = fu_in_s->fumate_p;
3655 
3656  NMG_GET_FU_NORMAL(normal, fu);
3657  if (reversed)
3658  VREVERSE(normal, normal);
3659 
3660  NMG_GET_FU_NORMAL(normal_in_s, fu_in_s);
3661 
3662  if (VDOT(normal, normal_in_s) < 0.0) {
3663  nmg_reverse_face_and_radials(fu_in_s, tol);
3664  }
3665  }
3666  }
3667 
3668  bu_ptbl_free(&reverse);
3669  bu_free((char *)trans_tbl, "translation table");
3670 
3671  nmg_km(tmp_m);
3672 }
3673 
3674 
3675 /**
3676  * This codes looks for situations as illustrated:
3677  *
3678  * *------->*-------->*--------->*
3679  * *<----------------------------*
3680  *
3681  * where one long edgeuse (the bottom one above) and two or more
3682  * shorter edgeuses (the top ones) are collinear and have the same
3683  * start and end vertices. The code breaks the longer edgeuse into
3684  * ones that can be radials of the shorter ones. Returns the number
3685  * of splits performed.
3686  */
3687 int
3688 nmg_break_long_edges(struct shell *s, const struct bn_tol *tol)
3689 {
3690  struct faceuse *fu;
3691  struct loopuse *lu;
3692  struct edgeuse *eu;
3693  int split_count=0;
3694 
3695  if (RTG.NMG_debug & DEBUG_BASIC)
3696  bu_log("nmg_break_long_edges(s = %p)\n", (void *)s);
3697 
3698  NMG_CK_SHELL(s);
3699  BN_CK_TOL(tol);
3700 
3701  /* look at every edgeuse in the shell */
3702  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
3703  NMG_CK_FACEUSE(fu);
3704 
3705  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
3706  NMG_CK_LOOPUSE(lu);
3707 
3708  /* skip loops of a single vertex */
3709  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
3710  continue;
3711 
3712  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
3713  struct vertexuse *vu;
3714 
3715  NMG_CK_EDGEUSE(eu);
3716 
3717  /* look for an edgeuse that terminates on this vertex */
3718  for (BU_LIST_FOR (vu, vertexuse, &eu->vu_p->v_p->vu_hd)) {
3719  struct edgeuse *eu1;
3720 
3721  NMG_CK_VERTEXUSE(vu);
3722 
3723  /* skip vertexuses that are not part of an edge */
3724  if (*vu->up.magic_p != NMG_EDGEUSE_MAGIC)
3725  continue;
3726 
3727  eu1 = vu->up.eu_p;
3728 
3729  /* don't consider the edge we already found! */
3730  if (eu1 == eu)
3731  continue;
3732 
3733  /* if it terminates at the same vertex as eu, skip it */
3734  if (eu1->eumate_p->vu_p->v_p == eu->eumate_p->vu_p->v_p)
3735  continue;
3736 
3737  /* get the mate (so it terminates at "vu") */
3738  eu1 = eu1->eumate_p;
3739 
3740  /* make sure it is collinear with "eu" */
3741  if (bn_3pts_collinear(eu->vu_p->v_p->vg_p->coord ,
3742  eu->eumate_p->vu_p->v_p->vg_p->coord ,
3743  eu1->vu_p->v_p->vg_p->coord, tol))
3744  {
3745  /* make sure we break the longer edge
3746  * and that the edges are in opposite directions */
3747  vect_t v0, v1;
3748 
3749  VSUB2(v0, eu->eumate_p->vu_p->v_p->vg_p->coord, eu->vu_p->v_p->vg_p->coord);
3750  VSUB2(v1, eu1->eumate_p->vu_p->v_p->vg_p->coord, eu1->vu_p->v_p->vg_p->coord);
3751 
3752  if (MAGSQ(v0) > MAGSQ(v1) && VDOT(v0, v1) < 0.0) {
3753  if (RTG.NMG_debug & DEBUG_BASIC) {
3754  bu_log("Breaking edge from (%f %f %f) to (%f %f %f) at (%f %f %f)\n",
3755  V3ARGS(eu->vu_p->v_p->vg_p->coord),
3756  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord),
3757  V3ARGS(eu1->vu_p->v_p->vg_p->coord));
3758  }
3759  (void) nmg_ebreak(eu1->vu_p->v_p, eu);
3760  split_count++;
3761  } else if (RTG.NMG_debug & DEBUG_BASIC) {
3762  bu_log("Not splitting collinear edges %p and %p:\n", (void *)eu, (void *)eu1);
3763  bu_log("\t(%f %f %f) -> (%f %f %f)\n",
3764  V3ARGS(eu->vu_p->v_p->vg_p->coord),
3765  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord));
3766  bu_log("\t(%f %f %f) -> (%f %f %f)\n",
3767  V3ARGS(eu1->vu_p->v_p->vg_p->coord),
3768  V3ARGS(eu1->eumate_p->vu_p->v_p->vg_p->coord));
3769  }
3770  }
3771  }
3772  }
3773  }
3774  }
3775  return split_count;
3776 }
3777 
3778 
3779 /**
3780  * Remove a loopuse from an existing face and construct a new face
3781  * from that loop
3782  *
3783  * Returns new faceuse as built by nmg_mf()
3784  */
3785 struct faceuse *
3787 {
3788  struct shell *s;
3789  struct faceuse *fu;
3790  struct loopuse *lu1;
3791  struct loopuse *lu_mate;
3792  int ot_same_loops=0;
3793 
3794  if (RTG.NMG_debug & DEBUG_BASIC)
3795  bu_log("nmg_mk_new_face_from_loop(lu = %p)\n", (void *)lu);
3796 
3797  NMG_CK_LOOPUSE(lu);
3798 
3799  if (*lu->up.magic_p != NMG_FACEUSE_MAGIC) {
3800  bu_log("nmg_mk_new_face_from_loop: loopuse is not in a faceuse\n");
3801  return (struct faceuse *)NULL;
3802  }
3803 
3804  fu = lu->up.fu_p;
3805  NMG_CK_FACEUSE(fu);
3806 
3807  s = fu->s_p;
3808  NMG_CK_SHELL(s);
3809 
3810  /* Count the number of exterior loops in this faceuse */
3811  for (BU_LIST_FOR (lu1, loopuse, &fu->lu_hd)) {
3812  NMG_CK_LOOPUSE(lu1);
3813  if (lu1->orientation == OT_SAME)
3814  ot_same_loops++;
3815  }
3816 
3817  if (ot_same_loops == 1 && lu->orientation == OT_SAME) {
3818  bu_log("nmg_mk_new_face_from_loop: cannot remove only exterior loop from faceuse\n");
3819  return (struct faceuse *)NULL;
3820  }
3821 
3822  lu_mate = lu->lumate_p;
3823 
3824  /* remove loopuse from faceuse */
3825  BU_LIST_DEQUEUE(&lu->l);
3826 
3827  /* remove its mate from faceuse mate */
3828  BU_LIST_DEQUEUE(&lu_mate->l);
3829 
3830  /* insert these loops in the shells list of wire loops
3831  * put the original loopuse at the head of the list
3832  * so that it will end up as the returned faceuse from "nmg_mf"
3833  */
3834  BU_LIST_INSERT(&s->lu_hd, &lu_mate->l);
3835  BU_LIST_INSERT(&s->lu_hd, &lu->l);
3836 
3837  /* set the "up" pointers to the shell */
3838  lu->up.s_p = s;
3839  lu_mate->up.s_p = s;
3840 
3841  /* Now make the new face */
3842  return nmg_mf(lu);
3843 }
3844 
3845 
3846 /* state for nmg_split_loops_into_faces */
3848 {
3849  long *flags; /* index based array of flags for model */
3850  const struct bn_tol *tol;
3851  int split; /* count of faces split */
3852 };
3853 
3854 
3855 void
3856 nmg_split_loops_handler(uint32_t *fu_p, void *sl_state, int UNUSED(unused))
3857 {
3858  struct faceuse *fu;
3859  struct nmg_split_loops_state *state;
3860  struct loopuse *lu;
3861  const struct bn_tol *tol;
3862  int otsame_loops=0;
3863  int otopp_loops=0;
3864 
3865  fu = (struct faceuse *)fu_p;
3866  NMG_CK_FACEUSE(fu);
3867 
3868  state = (struct nmg_split_loops_state *)sl_state;
3869  tol = state->tol;
3870 
3871  if (fu->orientation != OT_SAME)
3872  return;
3873 
3874  if (!NMG_INDEX_TEST_AND_SET(state->flags, fu)) return;
3875 
3876  NMG_INDEX_SET(state->flags, fu->fumate_p);
3877 
3878  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
3879  NMG_CK_LOOPUSE(lu);
3880 
3881  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
3882  continue;
3883 
3884  if (lu->orientation == OT_SAME)
3885  otsame_loops++;
3886  else if (lu->orientation == OT_OPPOSITE)
3887  otopp_loops++;
3888  else {
3889  bu_log("nmg_split_loops_into_faces: faceuse (%p) with %s loopuse (%p)\n",
3890  (void *)fu, nmg_orientation(lu->orientation), (void *)lu);
3891  return;
3892  }
3893  }
3894 
3895  /* if there is only one OT_SAME loop in this faceuse, nothing to do */
3896  if (otsame_loops == 1)
3897  return;
3898 
3899  if (otopp_loops && otsame_loops) {
3900  struct bu_ptbl inside_loops;
3901 
3902  bu_ptbl_init(&inside_loops, 64, " &inside_loops ");
3903 
3904  lu = BU_LIST_FIRST(loopuse, &fu->lu_hd);
3905  while (BU_LIST_NOT_HEAD(lu, &fu->lu_hd)) {
3906  struct faceuse *new_fu;
3907  struct loopuse *lu_next;
3908  struct loopuse *lu1;
3909  plane_t plane;
3910  int idx;
3911 
3912  lu_next = BU_LIST_PNEXT(loopuse, &lu->l);
3913 
3914  if (otsame_loops == 1) {
3915  /* done */
3916  bu_ptbl_free(&inside_loops);
3917  return;
3918  }
3919 
3920  if (lu->orientation != OT_SAME) {
3921  lu = lu_next;
3922  continue;
3923  }
3924 
3925  /* find OT_OPPOSITE loops for this exterior loop */
3926  for (BU_LIST_FOR (lu1, loopuse, &fu->lu_hd)) {
3927  struct loopuse *lu2;
3928  int hole_in_lu=1;
3929 
3930  if (lu1 == lu)
3931  continue;
3932 
3933  /* skip loops that are not OT_OPPOSITE */
3934  if (lu1->orientation != OT_OPPOSITE)
3935  continue;
3936 
3937  /* skip loops that are not within lu */
3938  if (nmg_classify_lu_lu(lu1, lu, tol) != NMG_CLASS_AinB)
3939  continue;
3940 
3941  /* lu1 is an OT_OPPOSITE loopuse within the OT_SAME lu
3942  * but lu1 may be within yet another loopuse that is
3943  * also within lu (nested loops)
3944  */
3945 
3946  /* look for another OT_SAME loopuse within lu */
3947  for (BU_LIST_FOR (lu2, loopuse, &fu->lu_hd)) {
3948 
3949  if (lu2 == lu || lu2 == lu1)
3950  continue;
3951 
3952  if (lu2->orientation != OT_SAME)
3953  continue;
3954 
3955  if (nmg_classify_lu_lu(lu2, lu, tol) == NMG_CLASS_AinB) {
3956  /* lu2 is within lu, does it contain lu1?? */
3957  if (nmg_classify_lu_lu(lu1, lu2, tol) == NMG_CLASS_AinB) {
3958  /* Yes, lu1 is within lu2, so lu1 is not
3959  * a hole in lu
3960  */
3961  hole_in_lu = 0;
3962  break;
3963  }
3964  }
3965  }
3966 
3967  if (hole_in_lu)
3968  bu_ptbl_ins(&inside_loops, (long *)lu1);
3969  }
3970 
3971  NMG_GET_FU_PLANE(plane, fu);
3972 
3973  new_fu = nmg_mk_new_face_from_loop(lu);
3974  if (new_fu) {
3975 
3976  nmg_face_g(new_fu, plane);
3977 
3978  for (idx=0; idx<BU_PTBL_END(&inside_loops); idx++) {
3979  lu1 = (struct loopuse *)BU_PTBL_GET(&inside_loops, idx);
3980  nmg_move_lu_between_fus(new_fu, fu, lu1);
3981  otopp_loops--;
3982  }
3983  nmg_face_bb(new_fu->f_p, tol);
3984  bu_ptbl_reset(&inside_loops);
3985  }
3986 
3987  otsame_loops--;
3988  lu = lu_next;
3989  }
3990  bu_ptbl_free(&inside_loops);
3991  } else if (otsame_loops) {
3992  /* all loops are of the same orientation, just make a face for every loop */
3993  int first=1;
3994  struct faceuse *new_fu;
3995 
3996  lu = BU_LIST_FIRST(loopuse, &fu->lu_hd);
3997  while (BU_LIST_NOT_HEAD(lu, &fu->lu_hd)) {
3998  struct loopuse *next_lu;
3999 
4000  next_lu = BU_LIST_PNEXT(loopuse, &lu->l);
4001 
4002  if (first)
4003  first = 0;
4004  else {
4005  plane_t plane;
4006 
4007  if (lu->orientation == OT_SAME) {
4008  NMG_GET_FU_PLANE(plane, fu);
4009  } else {
4010  NMG_GET_FU_PLANE(plane, fu->fumate_p);
4011  }
4012  new_fu = nmg_mk_new_face_from_loop(lu);
4013  if (new_fu) {
4014  nmg_face_g(new_fu, plane);
4015  nmg_face_bb(new_fu->f_p, tol);
4016  }
4017  }
4018 
4019  lu = next_lu;
4020  }
4021  } else {
4022  /* faceuse has only OT_OPPOSITE loopuses */
4023  bu_log("nmg_split_loops_into_faces: fu (%p) has only OT_OPPOSITE loopuses, ignored\n",
4024  (void *)fu);
4025  }
4026 }
4027 
4028 
4029 /**
4030  * Visits each faceuse and splits disjoint loops into separate faces.
4031  *
4032  * Returns the number of faces modified.
4033  */
4034 int
4035 nmg_split_loops_into_faces(uint32_t *magic_p, const struct bn_tol *tol)
4036 {
4037  struct model *m;
4038  struct nmg_split_loops_state sl_state;
4039  int count;
4040  static const struct nmg_visit_handlers htab = {NULL, NULL, NULL, NULL, NULL,
4041  NULL, NULL, NULL, NULL, nmg_split_loops_handler,
4042  NULL, NULL, NULL, NULL, NULL,
4043  NULL, NULL, NULL, NULL, NULL,
4044  NULL, NULL, NULL, NULL, NULL};
4045  /* htab.aft_faceuse = nmg_split_loops_handler; */
4046 
4047  if (RTG.NMG_debug & DEBUG_BASIC)
4048  bu_log("nmg_split_loops_into_faces(magic_p = %p)\n", (void *)magic_p);
4049 
4050  BN_CK_TOL(tol);
4051 
4052  m = nmg_find_model(magic_p);
4053  NMG_CK_MODEL(m);
4054 
4055  BN_CK_TOL(tol);
4056 
4057  sl_state.split = 0;
4058  sl_state.flags = (long *)bu_calloc(m->maxindex*2, sizeof(long), "nmg_split_loops_into_faces: flags");
4059  sl_state.tol = tol;
4060 
4061  nmg_visit(magic_p, &htab, (void *)&sl_state);
4062 
4063  count = sl_state.split;
4064 
4065  bu_free((char *)sl_state.flags, "nmg_split_loops_into_faces: flags");
4066 
4067  return count;
4068 }
4069 
4070 
4071 /**
4072  * Accepts one shell and breaks it to the minimum number of disjoint
4073  * shells.
4074  *
4075  * explicit returns:
4076  * # of resulting shells (1 indicates that nothing was done)
4077  *
4078  * implicit returns:
4079  * additional shells in the passed in shell's region.
4080  */
4081 int
4082 nmg_decompose_shell(struct shell *s, const struct bn_tol *tol)
4083 {
4084  int missed_faces;
4085  int no_of_shells=1;
4086  int shell_no=1;
4087  int i, j;
4088  struct model *m;
4089  struct nmgregion *r;
4090  struct shell *new_s;
4091  struct faceuse *fu;
4092  struct loopuse *lu;
4093  struct edgeuse *eu;
4094  struct edgeuse *eu1;
4095  struct faceuse *missed_fu = NULL;
4096  struct bu_ptbl stack;
4097  struct bu_ptbl shared_edges;
4098  long *flags;
4099 
4100  if (RTG.NMG_debug & DEBUG_BASIC)
4101  bu_log("nmg_decompose_shell(s = %p) START\n", (void *)s);
4102 
4103  NMG_CK_SHELL(s);
4104 
4105  BN_CK_TOL(tol);
4106 
4107  /* Make an index table to insure we visit each face once and only once */
4108  r = s->r_p;
4109  NMG_CK_REGION(r);
4110  m = r->m_p;
4111  NMG_CK_MODEL(m);
4112  flags = (long *)bu_calloc(m->maxindex*2, sizeof(long), "nmg_decompose_shell: flags");
4113 
4114  bu_ptbl_init(&stack, 64, " &stack ");
4115  bu_ptbl_init(&shared_edges, 64, " &shared_edges ");
4116 
4117  /* Need to be sure that every face has just one OT_SAME loop */
4118  (void)nmg_split_loops_into_faces(&s->l.magic, tol);
4119 
4120  /* get first faceuse from shell */
4121  fu = BU_LIST_FIRST(faceuse, &s->fu_hd);
4122  NMG_CK_FACEUSE(fu);
4123  if (fu->orientation != OT_SAME)
4124  fu = fu->fumate_p;
4125  if (fu->orientation != OT_SAME)
4126  bu_bomb("First face in shell has no OT_SAME uses!\n");
4127 
4128  /* put all edgeuses of first faceuse on the stack */
4129  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
4130  NMG_CK_LOOPUSE(lu);
4131  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
4132  continue;
4133 
4134  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
4135  /* build two lists, one of winged edges, the other not */
4136  if (nmg_radial_face_count(eu, s) > 2)
4137  bu_ptbl_ins_unique(&shared_edges, (long *)eu);
4138  else
4139  bu_ptbl_ins_unique(&stack, (long *)eu);
4140  }
4141  }
4142 
4143  /* Mark first faceuse and mate with shell number */
4144  NMG_INDEX_ASSIGN(flags, fu, shell_no);
4145  NMG_INDEX_ASSIGN(flags, fu->fumate_p, shell_no);
4146 
4147  /* now pop edgeuse of the stack and visit faces radial to edgeuse */
4148  while ((eu1 = nmg_pop_eu(&stack)) != (struct edgeuse *)NULL) {
4149  NMG_CK_EDGEUSE(eu1);
4150 
4151  /* move to the radial */
4152  eu = eu1->radial_p;
4153  NMG_CK_EDGEUSE(eu);
4154 
4155  /* make sure we stay within the intended shell */
4156  while (nmg_find_s_of_eu(eu) != s && eu != eu1 && eu != eu1->eumate_p)
4157  eu = eu->eumate_p->radial_p;
4158 
4159  if (eu == eu1 || eu == eu1->eumate_p)
4160  continue;
4161 
4162  fu = nmg_find_fu_of_eu(eu);
4163  NMG_CK_FACEUSE(fu);
4164 
4165  if (fu->orientation != OT_SAME)
4166  fu = fu->fumate_p;
4167 
4168  /* if this faceuse has already been visited, skip it */
4169  if (!NMG_INDEX_TEST(flags, fu)) {
4170  /* push all edgeuses of "fu" onto the stacks */
4171  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
4172  NMG_CK_LOOPUSE(lu);
4173 
4174  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
4175  continue;
4176 
4177  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
4178  /* build two lists, one of winged edges, the other not */
4179  if (nmg_radial_face_count(eu, s) > 2)
4180  bu_ptbl_ins_unique(&shared_edges, (long *)eu);
4181  else
4182  bu_ptbl_ins_unique(&stack, (long *)eu);
4183  }
4184  }
4185  /* Mark this faceuse and its mate with a shell number */
4186  NMG_INDEX_ASSIGN(flags, fu, shell_no);
4187  NMG_INDEX_ASSIGN(flags, fu->fumate_p, shell_no);
4188  }
4189  }
4190 
4191  /* count number of faces that were not visited */
4192  missed_faces = 0;
4193  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
4194  NMG_CK_FACEUSE(fu);
4195  if (fu->orientation == OT_SAME) {
4196  if (!NMG_INDEX_TEST(flags, fu)) {
4197  missed_faces++;
4198  missed_fu = fu;
4199  }
4200  }
4201  }
4202 
4203  if (!missed_faces) {
4204  /* nothing to do, just one shell */
4205  bu_free((char *)flags, "nmg_decompose_shell: flags ");
4206  bu_ptbl_free(&stack);
4207  bu_ptbl_free(&shared_edges);
4208  return no_of_shells;
4209  }
4210 
4211  while (missed_faces) {
4212  struct edgeuse *unassigned_eu = NULL;
4213  int *shells_at_edge;
4214  int new_shell_no=0;
4215 
4216  bu_ptbl_reset(&stack);
4217 
4218  /* Look at the list of shared edges to see if anything can be deduced */
4219  shells_at_edge = (int *)bu_calloc(no_of_shells+1, sizeof(int), "nmg_decompose_shell: shells_at_edge");
4220 
4221  for (i=0; i<BU_PTBL_END(&shared_edges); i++) {
4222  int faces_at_edge=0;
4223  int k;
4224 
4225  /* Construct a list of shells for this edge.
4226  * shells_at_edge[i] is the number of edgeuses of this
4227  * edge that have been assigned to shell number i.
4228  * shells_at_edge[0] is the number of uses of this edge
4229  * that have not been assigned to any shell yet
4230  */
4231  for (k=0; k<=no_of_shells; k++)
4232  shells_at_edge[k] = 0;
4233 
4234  unassigned_eu = NULL;
4235 
4236  eu = (struct edgeuse *)BU_PTBL_GET(&shared_edges, i);
4237  NMG_CK_EDGEUSE(eu);
4238 
4239  eu1 = eu;
4240  do {
4241  struct faceuse *fu_of_eu;
4242 
4243  fu_of_eu = nmg_find_fu_of_eu(eu1);
4244 
4245  faces_at_edge++;
4246  if (!NMG_INDEX_TEST(flags, fu_of_eu))
4247  unassigned_eu = eu1;
4248  shells_at_edge[ NMG_INDEX_GET(flags, fu_of_eu) ]++;
4249 
4250  eu1 = nmg_next_radial_eu(eu1, s, 0);
4251  }
4252  while (eu1 != eu);
4253 
4254  if (shells_at_edge[0] == 1 && unassigned_eu) {
4255  /* Only one edgeuse at this edge is unassigned, should be
4256  * able to determine which shell it belongs in */
4257 
4258  for (j=1; j<=no_of_shells; j++) {
4259  if (shells_at_edge[j] == 1) {
4260  /* unassigned edgeuse should belong to shell j */
4261  new_shell_no = j;
4262  break;
4263  }
4264  }
4265  } else if (shells_at_edge[0] == 0) {
4266  /* all uses of this edge have been assigned
4267  * it can be deleted from the table
4268  */
4269  bu_ptbl_rm(&shared_edges, (long *)eu);
4270  }
4271  if (new_shell_no)
4272  break;
4273  }
4274 
4275  if (!new_shell_no) {
4276  /* Couldn't find a definite edgeuse to start a new shell
4277  * so use radial parity to pick an edgeuse that should not be
4278  * part of the same shell as ones already assigned
4279  */
4280  for (i=0; i<BU_PTBL_END(&shared_edges); i++) {
4281  struct faceuse *fu_of_eu1;
4282  int found_missed_face=0;
4283 
4284  eu = (struct edgeuse *)BU_PTBL_GET(&shared_edges, i);
4285  NMG_CK_EDGEUSE(eu);
4286 
4287  eu1 = eu;
4288  do {
4289  /* look for unassigned edgeuses */
4290  fu_of_eu1 = nmg_find_fu_of_eu(eu1);
4291  NMG_CK_FACEUSE(fu_of_eu1);
4292  if (!NMG_INDEX_TEST(flags, fu_of_eu1)) {
4293  struct edgeuse *eu2;
4294  struct faceuse *fu_of_eu2;
4295 
4296  /* look for a neighboring edgeuse that
4297  * has been assigned
4298  */
4299  eu2 = nmg_prev_radial_eu(eu1, s, 0);
4300  fu_of_eu2 = nmg_find_fu_of_eu(eu2);
4301  NMG_CK_FACEUSE(fu_of_eu2);
4302  if (NMG_INDEX_TEST(flags, fu_of_eu2)) {
4303  /* eu2 has been assigned
4304  * compare orientation parity
4305  */
4306  if (fu_of_eu2->orientation ==
4307  fu_of_eu1->orientation) {
4308  /* These should not be in the same
4309  * shell, so start a new shell
4310  * at this faceuse
4311  */
4312  missed_fu = fu_of_eu1;
4313  found_missed_face = 1;
4314  }
4315  }
4316  if (found_missed_face)
4317  break;
4318 
4319  eu2 = nmg_next_radial_eu(eu1, s, 0);
4320  fu_of_eu2 = nmg_find_fu_of_eu(eu2);
4321  NMG_CK_FACEUSE(fu_of_eu2);
4322  if (NMG_INDEX_TEST(flags, fu_of_eu2)) {
4323  /* eu2 has been assigned
4324  * compare orientation parity
4325  */
4326  if (fu_of_eu2->orientation ==
4327  fu_of_eu1->orientation) {
4328  /* These should not be in the same
4329  * shell, so start a new shell
4330  * at this faceuse
4331  */
4332  missed_fu = fu_of_eu1;
4333  found_missed_face = 1;
4334  }
4335  }
4336 
4337  }
4338  if (found_missed_face)
4339  break;
4340  eu1 = nmg_next_radial_eu(eu1, s, 0);
4341  }
4342  while (eu1 != eu);
4343 
4344  if (found_missed_face)
4345  break;
4346  }
4347  }
4348 
4349  bu_free((char *)shells_at_edge, "nmg_decompose_shell: shells_at_edge");
4350 
4351  /* make a new shell number */
4352  if (new_shell_no) {
4353  shell_no = new_shell_no;
4354  fu = nmg_find_fu_of_eu(unassigned_eu);
4355  } else {
4356  shell_no = (++no_of_shells);
4357  NMG_CK_FACEUSE(missed_fu);
4358  fu = missed_fu;
4359  }
4360 
4361  if (fu->orientation != OT_SAME)
4362  fu = fu->fumate_p;
4363 
4364  if (!NMG_INDEX_TEST(flags, fu)) {
4365  /* move this missed face to the new shell */
4366 
4367  /* push all edgeuses of "fu" onto the stack */
4368  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
4369  NMG_CK_LOOPUSE(lu);
4370 
4371  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
4372  continue;
4373 
4374  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
4375  /* build two lists, one of winged edges, the other not */
4376  if (nmg_radial_face_count(eu, s) > 2)
4377  bu_ptbl_ins_unique(&shared_edges, (long *)eu);
4378  else
4379  bu_ptbl_ins_unique(&stack, (long *)eu);
4380  }
4381  }
4382 
4383  /* Mark this faceuse with a shell number */
4384  NMG_INDEX_ASSIGN(flags, fu, shell_no);
4385  NMG_INDEX_ASSIGN(flags, fu->fumate_p, shell_no);
4386 
4387  } else
4388  bu_bomb("nmg_decompose_shell: Missed face wasn't missed???\n");
4389 
4390  /* now pop edgeuse of the stack and visit faces radial to edgeuse */
4391  while ((eu1 = nmg_pop_eu(&stack)) != (struct edgeuse *)NULL) {
4392  NMG_CK_EDGEUSE(eu1);
4393 
4394  /* move to the radial */
4395  eu = eu1->radial_p;
4396  NMG_CK_EDGEUSE(eu);
4397 
4398  /* stay within the original shell */
4399  while (nmg_find_s_of_eu(eu) != s && eu != eu1 && eu != eu1->eumate_p)
4400  eu = eu->eumate_p->radial_p;
4401 
4402  if (eu == eu1 || eu == eu1->eumate_p)
4403  continue;
4404 
4405  fu = nmg_find_fu_of_eu(eu);
4406  NMG_CK_FACEUSE(fu);
4407 
4408  /* if this face has already been visited, skip it */
4409  if (!NMG_INDEX_TEST(flags, fu)) {
4410  /* push all edgeuses of "fu" onto the stack */
4411  for (BU_LIST_FOR (lu, loopuse, &fu->lu_hd)) {
4412  NMG_CK_LOOPUSE(lu);
4413  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
4414  continue;
4415  for (BU_LIST_FOR (eu, edgeuse, &lu->down_hd)) {
4416  /* build two lists, one of winged edges, the other not */
4417  if (nmg_radial_face_count(eu, s) > 2)
4418  bu_ptbl_ins_unique(&shared_edges, (long *)eu);
4419  else
4420  bu_ptbl_ins_unique(&stack, (long *)eu);
4421  }
4422  }
4423 
4424  /* Mark this faceuse with a shell number */
4425  NMG_INDEX_ASSIGN(flags, fu, shell_no);
4426  NMG_INDEX_ASSIGN(flags, fu->fumate_p, shell_no);
4427 
4428  }
4429  }
4430 
4431  /* count number of faces that were not visited */
4432  missed_faces = 0;
4433  for (BU_LIST_FOR (fu, faceuse, &s->fu_hd)) {
4434  NMG_CK_FACEUSE(fu);
4435  if (fu->orientation == OT_SAME) {
4436  if (!NMG_INDEX_TEST(flags, fu)) {
4437  missed_faces++;
4438  missed_fu = fu;
4439  }
4440  }
4441  }
4442  }
4443 
4444  /* Now build the new shells */
4445  for (shell_no=2; shell_no<=no_of_shells; shell_no++) {
4446  struct bu_ptbl faces;
4447 
4448  bu_ptbl_init(&faces, 64, "Faces ptbl for nmg_decompose_shell");
4449 
4450  /* Make a shell */
4451  new_s = nmg_msv(r);
4452  (void)nmg_kvu(new_s->vu_p);
4453 
4454  /* Move faces marked with this shell_no to this shell */
4455  fu = BU_LIST_FIRST(faceuse, &s->fu_hd);
4456  while (BU_LIST_NOT_HEAD(fu, &s->fu_hd)) {
4457  struct faceuse *next_fu;
4458 
4459  next_fu = BU_LIST_NEXT(faceuse, &fu->l);
4460  while (BU_LIST_NOT_HEAD(next_fu, &s->fu_hd) && next_fu->orientation != OT_SAME)
4461  next_fu = BU_LIST_NEXT(faceuse, &next_fu->l);
4462 
4463  if (fu->orientation != OT_SAME) {
4464  fu = next_fu;
4465  continue;
4466  }
4467 
4468  if (NMG_INDEX_GET(flags, fu) == shell_no) {
4469  nmg_mv_fu_between_shells(new_s, s, fu);
4470  bu_ptbl_ins(&faces, (long *)fu);
4471  }
4472 
4473  fu = next_fu;
4474  }
4475 
4476  nmg_gluefaces((struct faceuse **)BU_PTBL_BASEADDR(&faces), BU_PTBL_END(&faces), tol);
4477  bu_ptbl_free(&faces);
4478  nmg_shell_a(new_s, tol);
4479  }
4480  bu_free((char *)flags, "nmg_decompose_shell: flags ");
4481  bu_ptbl_free(&stack);
4482  bu_ptbl_free(&shared_edges);
4483 
4484  nmg_rebound(m , tol);
4485 
4486  if (RTG.NMG_debug & DEBUG_BASIC)
4487  bu_log("nmg_decompose_shell END (%d shells)\n", no_of_shells);
4488 
4489  return no_of_shells;
4490 }
4491 
4492 
4493 /**
4494  * Store an NMG model as a separate .g file, for later examination.
4495  * Don't free the model, as the caller may still have uses for it.
4496  *
4497  * NON-PARALLEL because of rt_uniresource.
4498  */
4499 void
4500 nmg_stash_model_to_file(const char *filename, const struct model *m, const char *title)
4501 {
4502  struct rt_wdb *fp;
4503  struct rt_db_internal intern;
4504  struct bu_external ext;
4505  int ret;
4506  int flags;
4507  char *name="error.s";
4508 
4509  bu_log("nmg_stash_model_to_file('%s', %p, %s)\n", filename, (void *)m, title);
4510 
4511  NMG_CK_MODEL(m);
4512  nmg_vmodel(m);
4513 
4514  if ((fp = wdb_fopen(filename)) == NULL) {
4515  perror(filename);
4516  return;
4517  }
4518 
4519  RT_DB_INTERNAL_INIT(&intern);
4520  intern.idb_major_type = DB5_MAJORTYPE_BRLCAD;
4521  intern.idb_type = ID_NMG;
4522  intern.idb_meth = &OBJ[ID_NMG];
4523  intern.idb_ptr = (void *)m;
4524 
4525  if (db_version(fp->dbip) < 5) {
4526  BU_EXTERNAL_INIT(&ext);
4527  ret = intern.idb_meth->ft_export4(&ext, &intern, 1.0, fp->dbip, &rt_uniresource);
4528  if (ret < 0) {
4529  bu_log("rt_db_put_internal(%s): solid export failure\n",
4530  name);
4531  ret = -1;
4532  goto out;
4533  }
4534  db_wrap_v4_external(&ext, name);
4535  } else {
4536  if (rt_db_cvt_to_external5(&ext, name, &intern, 1.0, fp->dbip, &rt_uniresource, intern.idb_major_type) < 0) {
4537  bu_log("wdb_export4(%s): solid export failure\n",
4538  name);
4539  ret = -2;
4540  goto out;
4541  }
4542  }
4543  BU_CK_EXTERNAL(&ext);
4544 
4545  flags = db_flags_internal(&intern);
4546  ret = wdb_export_external(fp, &ext, name, flags, intern.idb_type);
4547 out:
4548  bu_free_external(&ext);
4549  wdb_close(fp);
4550 
4551  bu_log("nmg_stash_model_to_file(): wrote error.s to '%s'\n",
4552  filename);
4553 }
4554 
4555 
4556 /* state for nmg_unbreak_edge */
4558 {
4559  long *flags; /* index based array of flags for model */
4560  int unbroken; /* count of edges mended */
4561 };
4562 
4563 
4564 /**
4565  * edgeuse visit routine for nmg_unbreak_region_edges.
4566  *
4567  * checks if edgeuse "eu" and its successor are candidates to be
4568  * unbroken. looks for two consecutive edgeuses sharing the same
4569  * edge geometry. Checks that the middle vertex has no other uses,
4570  * and, if so, kills the second edgeuse. Also moves the vu of the
4571  * first edgeuse mate to the vu of the killed edgeuse mate.
4572  */
4573 void
4574 nmg_unbreak_handler(uint32_t *eup, void *state, int UNUSED(unused))
4575 {
4576  struct edgeuse *eu1, *eu2;
4577  struct edge *e;
4578  struct edge_g_lseg *eg;
4579  struct nmg_unbreak_state *ub_state;
4580  struct vertex *vb;
4581  int ret;
4582 
4583  eu1 = (struct edgeuse *)eup;
4584  NMG_CK_EDGEUSE(eu1);
4585 
4586  ub_state = (struct nmg_unbreak_state *)state;
4587 
4588  /* there is a temptation to do a NMG_INDEX_SET(ub_state->flags, eu1->eumate_p)
4589  * here to avoid looking at this edgeuse's mate, but since we are only looking
4590  * forward, we must look at ALL edgeuses
4591  */
4592 
4593  /* make sure we only visit this edgeuse once */
4594  if (!NMG_INDEX_TEST_AND_SET(ub_state->flags, eu1)) return;
4595 
4596  e = eu1->e_p;
4597  NMG_CK_EDGE(e);
4598 
4599  eg = eu1->g.lseg_p;
4600  if (!eg) {
4601  bu_log("nmg_unbreak_handler: no geometry for edge %p\n", (void *)e);
4602  return;
4603  }
4604  NMG_CK_EDGE_G_EITHER(eg);
4605 
4606  /* if the edge geometry doesn't have at least two uses, this
4607  * is not a candidate for unbreaking */
4608  if (bu_list_len(&eg->eu_hd2) < 2*2) {
4609  /* bu_log("nmg_unbreak_handler: usage < 4\n"); */
4610  return;
4611  }
4612 
4613  /* Check for two consecutive uses, by looking forward. */
4614  eu2 = BU_LIST_PNEXT_CIRC(edgeuse, eu1);
4615  NMG_CK_EDGEUSE(eu2);
4616  if (eu2->g.lseg_p != eg) {
4617  /* Can't look backward here, or nmg_unbreak_edge()
4618  * will be asked to kill *this* edgeuse, which
4619  * will blow our caller's mind.
4620  */
4621  /* bu_log("nmg_unbreak_handler: eu1 edge geom not shared with eu2\n"); */
4622  return;
4623  }
4624  vb = eu2->vu_p->v_p;
4625  NMG_CK_VERTEX(vb);
4626 
4627  /* at this point, the situation is:
4628 
4629  eu1 eu2
4630  *----------->*----------->*
4631  A------------B------------C
4632  *<-----------*<-----------*
4633  eu1mate eu2mate
4634  */
4635  ret = nmg_unbreak_edge(eu1);
4636  if (ret != 0) return;
4637 
4638  /* keep a count of unbroken edges */
4639  ub_state->unbroken++;
4640 }
4641 
4642 
4643 /**
4644  * Uses the visit handler to call nmg_unbreak_handler for each edgeuse
4645  * below the region (or any other NMG element).
4646  *
4647  * returns the number of edges mended
4648  */
4649 int
4651 {
4652  struct model *m;
4653  struct nmg_unbreak_state ub_state;
4654  int count;
4655  static const struct nmg_visit_handlers htab = {NULL, NULL, NULL, NULL, NULL,
4656  NULL, NULL, NULL, NULL, NULL,
4657  NULL, NULL, NULL, NULL, NULL,
4658  NULL, NULL, nmg_unbreak_handler, NULL, NULL,
4659  NULL, NULL, NULL, NULL, NULL};
4660  /* htab.aft_edgeuse = nmg_unbreak_handler; */
4661 
4662  if (RTG.NMG_debug & DEBUG_BASIC)
4663  bu_log("nmg_unbreak_region_edges(magic_p = %p)\n", (void *)magic_p);
4664 
4665  m = nmg_find_model(magic_p);
4666  NMG_CK_MODEL(m);
4667 
4668  ub_state.unbroken = 0;
4669  ub_state.flags = (long *)bu_calloc(m->maxindex*2, sizeof(long), "nmg_unbreak_region_edges: flags");
4670 
4671  nmg_visit(magic_p, &htab, (void *)&ub_state);
4672 
4673  count = ub_state.unbroken;
4674 
4675  bu_free((char *)ub_state.flags, "nmg_unbreak_region_edges: flags");
4676 
4677  return count;
4678 }
4679 
4680 
4681 /**
4682  * Move a shell from one nmgregion to another. Will bomb if shell and
4683  * region aren't in the same model.
4684  *
4685  * returns:
4686  * 0 - all is well
4687  * 1 - nmgregion that gave up the shell is now empty!
4688  */
4689 int
4690 nmg_mv_shell_to_region(struct shell *s, struct nmgregion *r)
4691 {
4692  int ret_val;
4693 
4694  if (RTG.NMG_debug & DEBUG_BASIC)
4695  bu_log("nmg_mv_shell_to_region(s=%p, r=%p)\n", (void *)s, (void *)r);
4696 
4697  NMG_CK_SHELL(s);
4698  NMG_CK_REGION(r);
4699 
4700  if (s->r_p == r) {
4701  bu_log("nmg_mv_shell_to_region: Attempt to move shell to region it is already in\n");
4702  return 0;
4703  }
4704 
4705  if (nmg_find_model(&s->l.magic) != nmg_find_model(&r->l.magic))
4706  bu_bomb("nmg_mv_shell_to_region: Cannot move shell to a different model\n");
4707 
4708  BU_LIST_DEQUEUE(&s->l);
4709  if (BU_LIST_IS_EMPTY(&s->r_p->s_hd))
4710  ret_val = 1;
4711  else
4712  ret_val = 0;
4713 
4714  BU_LIST_APPEND(&r->s_hd, &s->l);
4715 
4716  s->r_p = r;
4717 
4718  return ret_val;
4719 }
4720 
4721 
4722 /**
4723  * Find all faces that contain vertex "new_v" Put them in a bu_ptbl
4724  * "faces"
4725  *
4726  * returns:
4727  * the number of faces that contain the vertex
4728  *
4729  * and fills in the table with the faces. Counts edges at this vertex
4730  * where radial is mate (free_edges)
4731  */
4732 int
4733 nmg_find_isect_faces(const struct vertex *new_v, struct bu_ptbl *faces, int *free_edges, const struct bn_tol *tol)
4734 {
4735  struct faceuse *fu;
4736  struct face_g_plane *fg;
4737  struct vertexuse *vu;
4738  int i;
4739  int unique;
4740 
4741  if (RTG.NMG_debug & DEBUG_BASIC)
4742  bu_log("nmg_find_isect_faces(new_v=%p, faces=%p)\n", (void *)new_v, (void *)faces);
4743 
4744  NMG_CK_VERTEX(new_v);
4745  BN_CK_TOL(tol);
4746  BU_CK_PTBL(faces);
4747 
4748  /* loop through vertex's vu list */
4749  for (BU_LIST_FOR (vu, vertexuse, &new_v->vu_hd)) {
4750  NMG_CK_VERTEXUSE(vu);
4751  fu = nmg_find_fu_of_vu(vu);
4752  if (fu->orientation != OT_SAME)
4753  continue;
4754 
4755  NMG_CK_FACEUSE(fu);
4756  fg = fu->f_p->g.plane_p;
4757 
4758  /* check if this face is different from the ones on list */
4759  unique = 1;
4760  for (i=0; i<BU_PTBL_END(faces); i++) {
4761  struct face *fp;
4762 
4763  fp = (struct face *)BU_PTBL_GET(faces, i);
4764  if (fp->g.plane_p == fg || bn_coplanar(fg->N, fp->g.plane_p->N, tol) > 0) {
4765  unique = 0;
4766  break;
4767  }
4768  }
4769 
4770  /* if it is not already on the list, add it */
4771  if (unique) {
4772  struct edgeuse *eu1;
4773 
4774  bu_ptbl_ins(faces, (long *)fu->f_p);
4775  /* Count the number of free edges containing new_v */
4776 
4777  if (*vu->up.magic_p != NMG_EDGEUSE_MAGIC)
4778  continue;
4779 
4780  eu1 = vu->up.eu_p;
4781  if (eu1->eumate_p == eu1->radial_p)
4782  (*free_edges)++;
4783  else {
4784  eu1 = BU_LIST_PPREV_CIRC(edgeuse, eu1);
4785  if (eu1->eumate_p == eu1->radial_p)
4786  (*free_edges)++;
4787  }
4788  }
4789  }
4790  return BU_PTBL_END(faces);
4791 }
4792 
4793 
4794 /**
4795  * given a vertex and a list of faces (not more than three) that
4796  * should intersect at the vertex, calculate a new location for the
4797  * vertex.
4798  *
4799  * returns:
4800  * 0 - if everything is OK
4801  * 1 - failure
4802  *
4803  * and modifies the geometry of the vertex to the new location
4804  */
4805 int
4806 nmg_simple_vertex_solve(struct vertex *new_v, const struct bu_ptbl *faces, const struct bn_tol *tol)
4807 {
4808  struct vertex_g *vg;
4809  int failed=0;
4810 
4811  if (RTG.NMG_debug & DEBUG_BASIC) {
4812  struct face *f;
4813  struct faceuse *fu;
4814  plane_t pl;
4815  int i;
4816 
4817  bu_log("nmg_simple_vertex_solve(new_v=%p, %ld faces)\n",
4818  (void *)new_v, BU_PTBL_END(faces));
4819 
4820  for (i = 0; i < BU_PTBL_END(faces); i++) {
4821  f = (struct face *)BU_PTBL_GET(faces, i);
4822  fu = f->fu_p;
4823  if (fu->orientation != OT_SAME)
4824  fu = fu->fumate_p;
4825  if (fu->orientation != OT_SAME)
4826  bu_log("\tface (%p) has no OT_SAME use\n", (void *)f);
4827  NMG_GET_FU_PLANE(pl, fu);
4828  bu_log("\t#%d: %g %g %g %g\n", i, V4ARGS(pl));
4829  }
4830  }
4831 
4832  NMG_CK_VERTEX(new_v);
4833  BU_CK_PTBL(faces);
4834  BN_CK_TOL(tol);
4835 
4836  vg = new_v->vg_p;
4837  NMG_CK_VERTEX_G(vg);
4838 
4839  switch (BU_PTBL_END(faces)) {
4840  struct face *fp1, *fp2, *fp3;
4841  plane_t pl1;
4842  fastf_t vert_move_len;
4843  fastf_t pl_dot;
4844 
4845  case 0:
4846  bu_log("nmg_simple_vertex_solve: vertex not in any face planes!\n");
4847  failed = 1;
4848  break;
4849 
4850  case 1: /* just move the vertex to the plane */
4851  fp1 = (struct face *)BU_PTBL_GET(faces, 0);
4852  vert_move_len = DIST_PT_PLANE(vg->coord, fp1->g.plane_p->N);
4853  VJOIN1(vg->coord, vg->coord, -vert_move_len, fp1->g.plane_p->N);
4854  break;
4855 
4856  case 2:
4857  fp1 = (struct face *)BU_PTBL_GET(faces, 0);
4858  fp2 = (struct face *)BU_PTBL_GET(faces, 1);
4859 
4860  pl_dot = VDOT(fp1->g.plane_p->N, fp2->g.plane_p->N);
4861  if (NEAR_EQUAL(pl_dot, 1.0, tol->perp) || NEAR_EQUAL(pl_dot, -1.0, tol->perp)) {
4862  vect_t move_vect;
4863 
4864  /* treat as a single plane */
4865  vert_move_len = (DIST_PT_PLANE(vg->coord, fp1->g.plane_p->N)
4866  + DIST_PT_PLANE(vg->coord, fp2->g.plane_p->N))/2.0;
4867  VADD2(move_vect, fp1->g.plane_p->N, fp2->g.plane_p->N);
4868  VUNITIZE(move_vect);
4869  VJOIN1(vg->coord, vg->coord, -vert_move_len, move_vect);
4870  } else {
4871  /* create a third plane perpendicular to first two */
4872 
4873  VCROSS(pl1, fp1->g.plane_p->N, fp2->g.plane_p->N);
4874 
4875  VUNITIZE(pl1);
4876  pl1[W] = VDOT(vg->coord, pl1);
4877  if (bn_mkpoint_3planes(vg->coord, fp1->g.plane_p->N, fp2->g.plane_p->N, pl1)) {
4878  bu_log("nmg_simple_vertex_solve: Cannot find new coords for two planes\n");
4879  bu_log("\tplanes are (%f %f %f %f) and (%f %f %f %f)\n",
4880  V4ARGS(fp1->g.plane_p->N) ,
4881  V4ARGS(fp2->g.plane_p->N));
4882  bu_log("\tcalculated third plane is (%f %f %f %f)\n", V4ARGS(pl1));
4883  failed = 1;
4884  break;
4885  }
4886  }
4887  break;
4888 
4889  case 3: /* just intersect the three planes */
4890  fp1 = (struct face *)BU_PTBL_GET(faces, 0);
4891  fp2 = (struct face *)BU_PTBL_GET(faces, 1);
4892  fp3 = (struct face *)BU_PTBL_GET(faces, 2);
4893  if (bn_mkpoint_3planes(vg->coord, fp1->g.plane_p->N, fp2->g.plane_p->N, fp3->g.plane_p->N)) {
4894  bu_log("nmg_simple_vertex_solve: failed for 3 planes:\n");
4895  bu_log("\t(%f %f %f %f)\n", V4ARGS(fp1->g.plane_p->N));
4896  bu_log("\t(%f %f %f %f)\n", V4ARGS(fp2->g.plane_p->N));
4897  bu_log("\t(%f %f %f %f)\n", V4ARGS(fp3->g.plane_p->N));
4898  failed = 1;
4899  break;
4900  }
4901  break;
4902  default:
4903  failed = 1;
4904  bu_log("nmg_simple_vertex_solve: Called for a complex vertex\n");
4905  break;
4906  }
4907 
4908  if (failed)
4909  bu_log("nmg_simple_vertex_solve: Failed to determine new coordinates for vertex at (%f %f %f)\n",
4910  V3ARGS(new_v->vg_p->coord));
4911  else if (RTG.NMG_debug & DEBUG_BASIC)
4912  bu_log("nmg_simple_vertex_solve: new coords = (%f %f %f)\n",
4913  V3ARGS(new_v->vg_p->coord));
4914 
4915  return failed;
4916 }
4917 
4918 
4919 /**
4920  * Check all uses of a vertex to see if it lies within tolerance of
4921  * all faces where it is used
4922  *
4923  * returns:
4924  * 0 - All is well
4925  * 1 - vertex is off face plane by at least tol->dist for at least one face
4926  */
4927 int
4928 nmg_ck_vert_on_fus(const struct vertex *v, const struct bn_tol *tol)
4929 {
4930  struct faceuse *fu;
4931  struct vertexuse *vu;
4932  fastf_t max_dist = 0.0;
4933  fastf_t dist = 0.0;
4934  int ret_val = 0;
4935  plane_t pl;
4936 
4937  NMG_CK_VERTEX(v);
4938  BN_CK_TOL(tol);
4939  NMG_CK_VERTEX_G(v->vg_p);
4940 
4941  for (BU_LIST_FOR (vu, vertexuse, &v->vu_hd)) {
4942  /* nmg_ck_vertexuse called within nmg_find_fu_of_vu,
4943  * so do not need to call here
4944  */
4945  fu = nmg_find_fu_of_vu(vu);
4946  if (!fu) {
4947  continue;
4948  }
4949  /* nmg_ck_faceuse, nmg_ck_face and nmg_ck_face_g_plane
4950  * are called within the nmg_get_fu_plane macro so
4951  * do not need to call them here
4952  */
4953  NMG_GET_FU_PLANE(pl, fu);
4954  dist = fabs(DIST_PT_PLANE(v->vg_p->coord, pl));
4955  if (dist > tol->dist) {
4956  ret_val = 1;
4957  if (dist > max_dist) {
4958  max_dist = dist;
4959  }
4960  bu_log("nmg_ck_vert_on_fus: v=%p vu=%p (%f %f %f) is %g from\n\tfaceuse %p f %p\n",
4961  (void *)v, (void *)vu, V3ARGS(v->vg_p->coord),
4962  dist,
4963  (void *)fu, (void *)fu->f_p);
4964  }
4965  }
4966 
4967  if (ret_val)
4968  bu_log("nmg_ck_vert_on_fus: v=%p (%f %f %f) max distance of %g from faceuses\n",
4969  (void *)v, V3ARGS(v->vg_p->coord), max_dist);
4970 
4971  return ret_val;
4972 }
4973 
4974 
4975 /* struct used by nmg_complex_vertex_solve
4976  * to store info about one edge
4977  * that contains the vertex in question
4978  */
4980 {
4981  struct faceuse *fu[2]; /* fu's that intersect at this edge */
4982  struct edgeuse *eu; /* edgeuse in fu[0] that emanates from vertex */
4983  point_t start; /* calculated start point of edge line */
4984  vect_t dir; /* calculated direction of edge line */
4985  point_t pt; /* a point on the edge a small distance from the vertex */
4986  int got_pt; /* flag indicating that the above point has been obtained */
4987  int free_edge; /* flag indicating that this is a free edge */
4988  struct vertex *vp; /* a vertex pointer for above point */
4989 };
4990 
4991 
4992 /**
4993  * debug printing of the table of intersect_fus structs used by
4994  * extruder
4995  */
4996 HIDDEN void
4997 nmg_pr_inter(const struct vertex *new_v, const struct bu_ptbl *int_faces)
4998 {
4999  int i;
5000  struct bn_tol tol;
5001 
5002  NMG_CK_VERTEX(new_v);
5003  BU_CK_PTBL(int_faces);
5004 
5005  tol.magic = BN_TOL_MAGIC;
5006  tol.dist = 0.0005;
5007  tol.dist_sq = tol.dist * tol.dist;
5008  tol.perp = 1e-6;
5009  tol.para = 1 - tol.perp;
5010 
5011  bu_log("\nint_faces at vertex %p (%f %f %f)\n", (void *)new_v, V3ARGS(new_v->vg_p->coord));
5012  for (i = 0; i < BU_PTBL_END(int_faces); i++) {
5013  struct intersect_fus *i_fus;
5014  struct face *fp1, *fp2;
5015  plane_t pl;
5016 
5017  i_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, i);
5018 
5019  bu_log("edge number %d, %p\n", i, (void *)i_fus);
5020  if (i_fus->fu[0])
5021  fp1 = i_fus->fu[0]->f_p;
5022  else
5023  fp1 = NULL;
5024  if (i_fus->fu[1]) {
5025  fp2 = i_fus->fu[1]->f_p;
5026  NMG_GET_FU_PLANE(pl, i_fus->fu[1]);
5027  } else
5028  fp2 = NULL;
5029 
5030  if (i_fus->fu[1])
5031  bu_log("\tfu1 = %p (face=%p), fu2 = %p (face=%p) (%f %f %f %f)\n",
5032  (void *)i_fus->fu[0], (void *)fp1, (void *)i_fus->fu[1], (void *)fp2,
5033  V4ARGS(pl));
5034  else
5035  bu_log("\tfu1 = %p (face=%p), fu2 = %p (face=%p)\n",
5036  (void *)i_fus->fu[0], (void *)fp1, (void *)i_fus->fu[1], (void *)fp2);
5037 
5038  if (i_fus->eu == NULL)
5039  bu_log("\teu = NULL\n");
5040  else if (i_fus->eu->vu_p == NULL)
5041  bu_log("\teu = %p, vu_p = NULL\n", (void *)i_fus->eu);
5042  else {
5043  struct faceuse *fu;
5044 
5045  bu_log("\teu = %p, from %p (%f %f %f) to %p (%f %f %f)\n",
5046  (void *)i_fus->eu,
5047  (void *)i_fus->eu->vu_p->v_p, V3ARGS(i_fus->eu->vu_p->v_p->vg_p->coord),
5048  (void *)i_fus->eu->eumate_p->vu_p->v_p, V3ARGS(i_fus->eu->eumate_p->vu_p->v_p->vg_p->coord));
5049  if (i_fus->fu[0]) {
5050  fu = nmg_find_fu_of_eu(i_fus->eu);
5051  if (fu != i_fus->fu[0])
5052  bu_log("****ERROR**** eu is not in fu1 it's in %p\n", (void *)fu);
5053  } else {
5054  fu = nmg_find_fu_of_eu(i_fus->eu);
5055  if (fu != i_fus->fu[1])
5056  bu_log("****ERROR**** eu is not in fu2, it's in %p\n", (void *)fu);
5057  }
5058  }
5059 
5060  bu_log("\tstart = (%f %f %f), dir = (%f %f %f)\n", V3ARGS(i_fus->start), V3ARGS(i_fus->dir));
5061  bu_log("\tpt = (%f %f %f)\n", V3ARGS(i_fus->pt));
5062  bu_log("\tfree_edge = %d\n", i_fus->free_edge);
5063  if (i_fus->eu && i_fus->eu->vu_p) {
5064  if (i_fus->eu->eumate_p != i_fus->eu->radial_p &&
5065  i_fus->free_edge)
5066  bu_log("****ERROR**** this is NOT a free edge\n");
5067  if (i_fus->eu->eumate_p == i_fus->eu->radial_p &&
5068  !i_fus->free_edge)
5069  bu_log("****ERROR**** This is a free edge\n");
5070  }
5071  if (i_fus->vp)
5072  bu_log("\tvp = %p (%f %f %f)\n", (void *)i_fus->vp, V3ARGS(i_fus->vp->vg_p->coord));
5073  else
5074  bu_log("\tvp = NULL\n");
5075  }
5076 }
5077 
5078 
5079 /**
5080  * Fill in the intersect_fus structures for edges around new_v. Does
5081  * not fill in "pt" or "vp".
5082  *
5083  * returns:
5084  * 0 - All is well
5085  * 1 - Failure
5086  */
5087 
5088 HIDDEN int
5089 nmg_get_edge_lines(struct vertex *new_v, struct bu_ptbl *int_faces, const struct bn_tol *tol)
5090 {
5091  struct vertex_g *vg;
5092  struct vertexuse *vu;
5093  struct edgeuse *eu, *eu1;
5094  struct faceuse *fu;
5095  struct model *m;
5096  struct nmgregion *r;
5097  struct bn_tol tol_tmp;
5098  int done = 0;
5099  int edge_no;
5100 
5101  NMG_CK_VERTEX(new_v);
5102  vg = new_v->vg_p;
5103  NMG_CK_VERTEX_G(vg);
5104  BN_CK_TOL(tol);
5105  BU_CK_PTBL(int_faces);
5106 
5107  if (RTG.NMG_debug & DEBUG_BASIC)
5108  bu_log("nmg_get_edge_lines(new_v=%p, int_faces=%p)\n", (void *)new_v, (void *)int_faces);
5109 
5110  /* A temporary tolerance struct for times when I don't want tolerancing */
5111  tol_tmp.magic = BN_TOL_MAGIC;
5112  tol_tmp.dist = 0.0;
5113  tol_tmp.dist_sq = 0.0;
5114  tol_tmp.perp = 0.0;
5115  tol_tmp.para = 1.0;
5116 
5117  m = nmg_find_model(&new_v->magic);
5118  NMG_CK_MODEL(m);
5119  r = BU_LIST_FIRST(nmgregion, &m->r_hd);
5120  NMG_CK_REGION(r);
5121  NMG_CK_REGION_A(r->ra_p);
5122 
5123  /* look for a dangling edge emanating from this vertex */
5124  eu1 = (struct edgeuse *)NULL;
5125  for (BU_LIST_FOR (vu, vertexuse, &new_v->vu_hd)) {
5126  if (*vu->up.magic_p != NMG_EDGEUSE_MAGIC)
5127  continue;
5128 
5129  eu = vu->up.eu_p->eumate_p;
5130  fu = nmg_find_fu_of_eu(eu);
5131  if (!fu)
5132  continue;
5133 
5134  if (fu->orientation != OT_SAME)
5135  continue;
5136 
5137  if (eu->eumate_p == eu->radial_p) {
5138  /* found a dangling edge, start processing here */
5139  plane_t pl;
5140  struct intersect_fus *i_fus;
5141 
5142  /* create and initialize an intersect_fus struct for this edge */
5143  BU_ALLOC(i_fus, struct intersect_fus);
5144  i_fus->fu[0] = NULL;
5145  i_fus->fu[1] = fu;
5146  i_fus->eu = eu;
5147  i_fus->vp = (struct vertex *)NULL;
5148  VSET(i_fus->pt, 0.0, 0.0, 0.0);
5149  i_fus->got_pt = 0;
5150  i_fus->free_edge = 1;
5151  eu1 = BU_LIST_PNEXT_CIRC(edgeuse, &eu->l);
5152 
5153  VSUB2(i_fus->dir, eu->vu_p->v_p->vg_p->coord, eu->eumate_p->vu_p->v_p->vg_p->coord);
5154  VUNITIZE(i_fus->dir);
5155  NMG_GET_FU_PLANE(pl, fu);
5156  VJOIN1(i_fus->start, vg->coord, (-DIST_PT_PLANE(vg->coord, pl)), pl);
5157 
5158  /* Save this info in the int_faces table */
5159  bu_ptbl_ins(int_faces, (long *)i_fus);
5160 
5161  break;
5162  }
5163  }
5164 
5165  if (!eu1) {
5166  int found_start=0;
5167 
5168  /* get the an edgeuse emanating from new_v */
5169  for (BU_LIST_FOR (vu, vertexuse, &new_v->vu_hd)) {
5170  NMG_CK_VERTEXUSE(vu);
5171  if (*vu->up.magic_p != NMG_EDGEUSE_MAGIC)
5172  continue;
5173 
5174  eu1 = vu->up.eu_p;
5175 
5176  fu = nmg_find_fu_of_eu(eu1);
5177  NMG_CK_FACEUSE(fu);
5178 
5179  if (fu->orientation == OT_SAME) {
5180  found_start = 1;
5181  break;
5182  }
5183  }
5184  if (!found_start) {
5185  bu_log("Cannot find edgeuse in OT_SAME faceuse starting at (%f %f %f)\n",
5186  V3ARGS(new_v->vg_p->coord));
5187  return 1;
5188  }
5189  }
5190 
5191  eu = eu1;
5192 
5193  /* loop through all the edges emanating from new_v */
5194  while (!done) {
5195  fastf_t dist;
5196  point_t start;
5197  vect_t dir;
5198  vect_t eu_dir;
5199  int ret_val;
5200  struct intersect_fus *i_fus;
5201  struct faceuse *fu1, *fu2;
5202 
5203  NMG_CK_EDGEUSE(eu);
5204 
5205  if (eu->vu_p->v_p != new_v) {
5206  /* This can happen if the faces of the shell are not properly
5207  * oriented such as might happen when an object is incorrectly
5208  * modelled in FASTGEN and run through the patch-g converter
5209  */
5210  bu_log("nmg_get_edge_lines: Bad solid!\n");
5211  for (edge_no=0; edge_no<BU_PTBL_END(int_faces); edge_no++) {
5212  i_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, edge_no);
5213 
5214  bu_free((char *)i_fus, "nmg_get_edge_lines: i_fus");
5215  }
5216  return 1;
5217  }
5218 
5219  /* get the direction of the original edge (away from the vertex) */
5220  VSUB2(eu_dir, eu->eumate_p->vu_p->v_p->vg_p->coord, eu->vu_p->v_p->vg_p->coord);
5221 
5222  /* get the two faces that intersect along this edge */
5223  fu1 = nmg_find_fu_of_eu(eu);
5224  fu2 = nmg_find_fu_of_eu(eu->radial_p);
5225 
5226  /* initialize the intersect structure for this edge */
5227  BU_ALLOC(i_fus, struct intersect_fus);
5228  i_fus->fu[0] = fu1;
5229  if (eu->radial_p == eu->eumate_p) {
5230  i_fus->fu[1] = (struct faceuse *)NULL;
5231  i_fus->free_edge = 1;
5232  done = 1;
5233  } else {
5234  i_fus->fu[1] = fu2;
5235  i_fus->free_edge = 0;
5236  }
5237  i_fus->eu = eu;
5238  i_fus->vp = (struct vertex *)NULL;
5239  VSET(i_fus->pt, 0.0, 0.0, 0.0);
5240  i_fus->got_pt = 0;
5241  VSET(i_fus->start, 0.0, 0.0, 0.0);
5242  VSET(i_fus->dir, 0.0, 0.0, 0.0);
5243 
5244  /* if edge is between loops of same face, don't calculate an edge line */
5245  if (fu1->f_p != fu2->f_p) {
5246  /* find the new edge line at the intersection of these two faces
5247  * the line is defined by start and dir */
5248 
5249  ret_val = bn_isect_2planes(start, dir,
5250  fu1->f_p->g.plane_p->N,
5251  fu2->f_p->g.plane_p->N,
5252  new_v->vg_p->coord,
5253  &tol_tmp);
5254  if (ret_val) {
5255  /* Cannot find line for this edge */
5256  bu_log("nmg_inside_vert: Cannot find new edge between two planes\n");
5257  bu_log("return from bn_isect_2planes is %d\n", ret_val);
5258  bu_log("\tplanes are (%f %f %f %f) and (%f %f %f %f)\n" ,
5259  V4ARGS(fu1->f_p->g.plane_p->N),
5260  V4ARGS(fu2->f_p->g.plane_p->N));
5261  bu_log("\tfus %p and %p, faces %p and %p\n" ,
5262  (void *)fu1, (void *)fu2, (void *)fu1->f_p, (void *)fu2->f_p);
5263  nmg_pr_fu_briefly(fu1, "fu1: ");
5264  nmg_pr_fu_briefly(fu2, "fu2: ");
5265  bu_bomb("Can't find plane intersection\n");
5266  }
5267  /* Make the start point at closest approach to old vertex */
5268  (void)bn_dist_pt3_line3(&dist, start, start, dir, new_v->vg_p->coord, tol);
5269 
5270  /* Make sure the calculated direction is away from the vertex */
5271  if (VDOT(eu_dir, dir) < 0.0)
5272  VREVERSE(dir, dir);
5273  VMOVE(i_fus->start, start);
5274  VMOVE(i_fus->dir, dir);
5275  } else if (i_fus->free_edge) {
5276  plane_t pl;
5277 
5278  /* for a dangling edge, use the same direction as the original edge
5279  * just move the start point to the new plane
5280  */
5281 
5282  NMG_GET_FU_PLANE(pl, fu1);
5283 
5284  VMOVE(i_fus->dir, eu_dir);
5285  VUNITIZE(i_fus->dir);
5286 
5287  VJOIN1(i_fus->start, vg->coord, (-DIST_PT_PLANE(vg->coord, pl)), pl);
5288 
5289  }
5290 
5291  /* Save this info in the int_faces table */
5292  bu_ptbl_ins(int_faces, (long *)i_fus);
5293 
5294  if (!done) {
5295  /* move on to the next edge emanating from new_v */
5296  eu = eu->radial_p;
5297  eu = BU_LIST_PNEXT_CIRC(edgeuse, eu);
5298  if (eu == eu1)
5299  done = 1;
5300  }
5301  }
5302  if (RTG.NMG_debug & DEBUG_BASIC) {
5303  bu_log("After getting edge lines:\n");
5304  nmg_pr_inter(new_v, int_faces);
5305  }
5306 
5307  return 0;
5308 }
5309 
5310 
5311 /**
5312  * Fill in the "pt" portion of the "intersect_fus" structure for edges
5313  * around new_v by calculating the intersection with neighboring edges
5314  * and selecting the furthest one from new_v.
5315  */
5316 HIDDEN int
5317 nmg_get_max_edge_inters(const struct vertex *new_v, struct bu_ptbl *int_faces, const struct bu_ptbl *faces, const struct bn_tol *tol)
5318 {
5319  struct model *m;
5320  struct nmgregion *r;
5321  int edge_no;
5322 
5323  if (RTG.NMG_debug & DEBUG_BASIC)
5324  bu_log("nmg_get_max_edge_inters(new_v = %p, %ld intersect_fus structs, %ld faces)\n",
5325  (void *)new_v, BU_PTBL_END(int_faces), BU_PTBL_END(faces));
5326 
5327  NMG_CK_VERTEX(new_v);
5328  BN_CK_TOL(tol);
5329  BU_CK_PTBL(int_faces);
5330 
5331  m = nmg_find_model(&new_v->magic);
5332  NMG_CK_MODEL(m);
5333  r = BU_LIST_FIRST(nmgregion, &m->r_hd);
5334  NMG_CK_REGION(r);
5335  NMG_CK_REGION_A(r->ra_p);
5336 
5337  /* loop through edges departing from new_v */
5338  for (edge_no=0; edge_no<BU_PTBL_END(int_faces); edge_no++) {
5339  struct intersect_fus *edge_fus, *other_fus;
5340  fastf_t max_dist, dist[2];
5341  int next_edge_no, prev_edge_no;
5342  int other_index;
5343 
5344  edge_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, edge_no);
5345 
5346  /* don't calculate intersect point for edge between two loops of same face */
5347  if (edge_fus->fu[0] && edge_fus->fu[1] &&
5348  edge_fus->fu[0]->f_p == edge_fus->fu[1]->f_p)
5349  continue;
5350 
5351  /* Find intersections with neighboring edges and keep the one
5352  * furthest up the edge
5353  */
5354  max_dist = (-MAX_FASTF);
5355 
5356  /* start with next edge */
5357  next_edge_no = edge_no + 1;
5358  if (next_edge_no == BU_PTBL_END(int_faces))
5359  next_edge_no = 0;
5360 
5361  other_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, next_edge_no);
5362 
5363  /* skip over edges between loops of same face */
5364  while (other_fus->fu[0] == other_fus->fu[1] && other_fus != edge_fus) {
5365  next_edge_no++;
5366  if (next_edge_no == BU_PTBL_END(int_faces))
5367  next_edge_no = 0;
5368 
5369  other_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, next_edge_no);
5370 
5371  }
5372 
5373  /* if we found another edge, calculate its intersection with the edge */
5374  if (other_fus != edge_fus) {
5375  if (!bn_dist_line3_line3(dist, edge_fus->start, edge_fus->dir, other_fus->start, other_fus->dir, tol)) {
5376  if (RTG.NMG_debug & DEBUG_BASIC)
5377  bu_log("Edge #%d intersects edge #%d at dist = %f\n", edge_no, next_edge_no, dist[0]);
5378  if (NEAR_ZERO(dist[0], tol->dist))
5379  dist[0] = 0.0;
5380  if (dist[0] > max_dist)
5381  max_dist = dist[0];
5382  }
5383  }
5384 
5385  /* now check the previous neighboring edge */
5386  prev_edge_no = edge_no - 1;
5387  if (prev_edge_no < 0)
5388  prev_edge_no = BU_PTBL_END(int_faces) - 1;
5389 
5390  other_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, prev_edge_no);
5391 
5392  while (other_fus->fu[0] == other_fus->fu[1] && other_fus != edge_fus) {
5393  prev_edge_no--;
5394  if (prev_edge_no < 0)
5395  prev_edge_no = BU_PTBL_END(int_faces) - 1;
5396 
5397  other_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, prev_edge_no);
5398  }
5399 
5400  if (other_fus != edge_fus) {
5401  if (bn_dist_line3_line3(dist, edge_fus->start, edge_fus->dir, other_fus->start, other_fus->dir, tol) >= 0) {
5402  if (RTG.NMG_debug & DEBUG_BASIC)
5403  bu_log("Edge #%d intersects edge #%d at dist = %f\n", edge_no, prev_edge_no, dist[0]);
5404  if (NEAR_ZERO(dist[0], tol->dist))
5405  dist[0] = 0.0;
5406  if (dist[0] > max_dist)
5407  max_dist = dist[0];
5408  }
5409  }
5410 
5411  if (max_dist < 0.0) {
5412  /* Now check for intersections with other planes */
5413  for (other_index=0; other_index<BU_PTBL_END(int_faces); other_index ++) {
5414  struct face *f;
5415 
5416  if (other_index == edge_no)
5417  continue;
5418 
5419  other_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, other_index);
5420 
5421  if (!other_fus->fu[0])
5422  continue;
5423 
5424  NMG_CK_FACEUSE(other_fus->fu[0]);
5425  f = other_fus->fu[0]->f_p;
5426 
5427  if (edge_fus->fu[0] && f == edge_fus->fu[0]->f_p)
5428  continue;
5429 
5430  if (edge_fus->fu[1] && f == edge_fus->fu[1]->f_p)
5431  continue;
5432 
5433  /* Do not intersect with a plane that this edge is parallel to */
5434  if (NEAR_ZERO(VDOT(f->g.plane_p->N, edge_fus->dir), tol->perp))
5435  continue;
5436 
5437  if (bn_isect_line3_plane(&dist[0], edge_fus->start, edge_fus->dir, f->g.plane_p->N, tol) > 1)
5438  continue;
5439 
5440  if (RTG.NMG_debug & DEBUG_BASIC)
5441  bu_log("Edge #%d intersects fu[0] from edge #%d at dist = %f\n", edge_no, other_index, dist[0]);
5442 
5443  if (NEAR_ZERO(dist[0], tol->dist))
5444  dist[0] = 0.0;
5445 
5446  if (dist[0] > max_dist)
5447  max_dist = dist[0];
5448  }
5449  }
5450 
5451  /* if any intersections have been found, save the point in edge_fus->pt */
5452  if (max_dist > (-MAX_FASTF)) {
5453  VJOIN1(edge_fus->pt, edge_fus->start, max_dist, edge_fus->dir);
5454  edge_fus->got_pt = 1;
5455  }
5456  }
5457 
5458  /* if no intersection was found, just use the edge-line start point */
5459  for (edge_no=0; edge_no < BU_PTBL_END(int_faces); edge_no++) {
5460  struct intersect_fus *edge_fus;
5461 
5462  edge_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, edge_no);
5463  if (!edge_fus->got_pt) {
5464  VMOVE(edge_fus->pt, edge_fus->start);
5465  }
5466  }
5467 
5468  if (RTG.NMG_debug & DEBUG_BASIC) {
5469  bu_log("After nmg_get_max_edge_inters:\n");
5470  nmg_pr_inter(new_v, int_faces);
5471  }
5472 
5473  return 0;
5474 }
5475 
5476 
5477 /**
5478  * eliminate "j_fus" from the table "int_faces" and adjust the info in
5479  * "i_fus". This is done when the "vp" vertices of the two structures
5480  * have been joined.
5481  */
5482 HIDDEN void
5483 nmg_fuse_inters(struct intersect_fus *i_fus, struct intersect_fus *j_fus, struct bu_ptbl *int_faces, const struct bn_tol *tol)
5484 {
5485  struct edgeuse *radial_eu;
5486  struct edgeuse *prev_eu;
5487 
5488  BU_CK_PTBL(int_faces);
5489  BN_CK_TOL(tol);
5490 
5491  if (RTG.NMG_debug & DEBUG_BASIC)
5492  bu_log("nmg_fuse_inters: i_fus=%p, j_fus=%p, %ld edges\n",
5493  (void *)i_fus, (void *)j_fus, BU_PTBL_END(int_faces));
5494 
5495  /* remember the radial edge of the structure to be deleted */
5496  radial_eu = j_fus->eu->radial_p;
5497 
5498  /* if the vertices have been joined prev_eu and j_fus->eu should be adjacent */
5499  prev_eu = BU_LIST_PPREV_CIRC(edgeuse, &j_fus->eu->l);
5500 
5501  if (EDGESADJ(prev_eu, j_fus->eu)) {
5502  nmg_keu(prev_eu);
5503  nmg_keu(j_fus->eu);
5504  } else
5505  bu_log("nmg_fuse_inter_verts: ERROR: can't find adjacent edges to kill\n");
5506 
5507  /* the other face for this edge is now j_fus->fu[1] */
5508  i_fus->fu[1] = j_fus->fu[1];
5509 
5510  /* if there are other faces along the edges that have been brought together
5511  * do a radial join
5512  */
5513  if (i_fus->fu[0] && j_fus->fu[1]) {
5514  if (RTG.NMG_debug & DEBUG_BASIC) {
5515  bu_log("radial join of eu's %p and %p\n", (void *)i_fus->eu, (void *)radial_eu);
5516  bu_log("\t%p to %p and %p to %p\n" ,
5517  (void *)i_fus->eu->vu_p->v_p, (void *)i_fus->eu->eumate_p->vu_p->v_p,
5518  (void *)radial_eu->vu_p->v_p, (void *)radial_eu->eumate_p->vu_p->v_p);
5519  }
5520  nmg_radial_join_eu(i_fus->eu, radial_eu, tol);
5521  }
5522 
5523  /* If this is a dangling edge, need to adjust the eu pointer */
5524  if (!i_fus->fu[0])
5525  i_fus->eu = radial_eu;
5526  NMG_CK_EDGEUSE(i_fus->eu);
5527 
5528  /* if the deleted structure was for a dangling edge,
5529  * then this edge is now dangling
5530  */
5531  if (j_fus->free_edge)
5532  i_fus->free_edge = 1;
5533 
5534  bu_ptbl_rm(int_faces, (long *)j_fus);
5535  bu_free((char *)j_fus, "nmg_split_edges_at_pts: j_fus ");
5536 
5537 }
5538 
5539 
5540 /**
5541  * Using the info in the table of intersect_fus structs, split the
5542  * edgeuse (eu) in each struct at the point (pt) store the new
5543  * vertices in the structure (vp) and assign the geometry.
5544  */
5545 HIDDEN void
5546 nmg_split_edges_at_pts(const struct vertex *new_v, struct bu_ptbl *int_faces, const struct bn_tol *tol)
5547 {
5548  int edge_no;
5549 
5550  if (RTG.NMG_debug & DEBUG_BASIC)
5551  bu_log("nmg_split_edges_at_pts(new_v = %p, %ld intersect_fus structs)\n",
5552  (void *)new_v, BU_PTBL_END(int_faces));
5553 
5554  BN_CK_TOL(tol);
5555  BU_CK_PTBL(int_faces);
5556  NMG_CK_VERTEX(new_v);
5557 
5558  /* loop through all edges departing from new_v */
5559  for (edge_no = 0; edge_no < BU_PTBL_END(int_faces); edge_no++) {
5560  struct intersect_fus *i_fus;
5561  struct edgeuse *new_eu;
5562 
5563  i_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, edge_no);
5564  if (i_fus == NULL)
5565  continue;
5566 
5567  /* skip edges between two loops of same face, for now */
5568  if (i_fus->fu[0] && i_fus->fu[1] && i_fus->fu[0]->f_p == i_fus->fu[1]->f_p)
5569  continue;
5570 
5571  if (bn_pt3_pt3_equal(new_v->vg_p->coord, i_fus->pt, tol)) {
5572  /* if pt is within tolerance of new_v, don't split the edge */
5573  i_fus->vp = (struct vertex *)NULL;
5574  VMOVE(i_fus->pt, new_v->vg_p->coord);
5575  VMOVE(i_fus->start, new_v->vg_p->coord);
5576  VSUB2(i_fus->dir, i_fus->eu->eumate_p->vu_p->v_p->vg_p->coord, i_fus->eu->vu_p->v_p->vg_p->coord);
5577  VUNITIZE(i_fus->dir);
5578  continue;
5579  }
5580  new_eu = nmg_esplit(i_fus->vp, i_fus->eu, 0);
5581  i_fus->vp = new_eu->vu_p->v_p;
5582 
5583  /* Need to keep track of correct eu in this case */
5584  if (i_fus->free_edge && !i_fus->fu[0])
5585  i_fus->eu = new_eu;
5586 
5587  /* Assign geometry to the new vertex */
5588  if (i_fus && !i_fus->vp->vg_p)
5589  nmg_vertex_gv(i_fus->vp, i_fus->pt);
5590  }
5591  if (RTG.NMG_debug & DEBUG_BASIC) {
5592  bu_log("After splitting edges:\n");
5593  nmg_pr_inter(new_v, int_faces);
5594  }
5595 
5596  /* Now take care of edges between two loops of same face */
5597  edge_no = 0;
5598  while (edge_no < BU_PTBL_END(int_faces)) {
5599  int next_edge_no;
5600  struct intersect_fus *i_fus, *j_fus;
5601 
5602  next_edge_no = edge_no + 1;
5603  if (next_edge_no == BU_PTBL_END(int_faces))
5604  next_edge_no = 0;
5605 
5606  i_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, edge_no);
5607  j_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, next_edge_no);
5608 
5609  /* look at all edges in the same face as i_fus->fu[1] */
5610  while (j_fus->fu[0] && j_fus->fu[1] &&
5611  j_fus->fu[0]->f_p == j_fus->fu[1]->f_p &&
5612  j_fus != i_fus)
5613  {
5614  /* if both edges are dangling, there is nothing to do */
5615  if (i_fus->free_edge && j_fus->free_edge)
5616  break;
5617 
5618  /* if we haven't assigned a vertex, skip this edge */
5619  if (!i_fus->vp)
5620  break;
5621 
5622  /* split the neighbor at the first structure's "vp"
5623  * this moves the neighboring edge's endpoint to
5624  * fall on the first edge.
5625  */
5626  (void) nmg_esplit(i_fus->vp, j_fus->eu, 0);
5627 
5628  /* now we can ignore this edge */
5629  nmg_fuse_inters(i_fus, j_fus, int_faces, tol);
5630 
5631  /* go to the next edge */
5632  if (next_edge_no == BU_PTBL_END(int_faces))
5633  next_edge_no = 0;
5634 
5635  j_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, next_edge_no);
5636 
5637  }
5638  edge_no++;
5639  }
5640  if (RTG.NMG_debug & DEBUG_BASIC) {
5641  bu_log("After loops of same face\n");
5642  nmg_pr_inter(new_v, int_faces);
5643  }
5644 }
5645 
5646 
5647 /**
5648  * kill all zero length edgeuses in faces around new_v
5649  */
5650 HIDDEN void
5651 nmg_remove_short_eus_inter(struct vertex *new_v, struct bu_ptbl *int_faces, const struct bn_tol *tol)
5652 {
5653  int edge_no;
5654  struct vertexuse *vu;
5655 
5656  NMG_CK_VERTEX(new_v);
5657  BU_CK_PTBL(int_faces);
5658  BN_CK_TOL(tol);
5659 
5660  if (RTG.NMG_debug & DEBUG_BASIC)
5661  bu_log("nmg_remove_short_eus: new_v=%p (%f %f %f), %ld edges\n",
5662  (void *)new_v, V3ARGS(new_v->vg_p->coord), BU_PTBL_END(int_faces));
5663 
5664  /* first join any of the "vp" in the intersect_fus structs that are
5665  * within tolerance of new-v
5666  */
5667  for (edge_no=0; edge_no<BU_PTBL_END(int_faces); edge_no++) {
5668  struct intersect_fus *edge_fus;
5669 
5670  edge_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, edge_no);
5671 
5672  if (!edge_fus->vp)
5673  continue;
5674 
5675  if (!bn_pt3_pt3_equal(new_v->vg_p->coord, edge_fus->vp->vg_p->coord, tol))
5676  continue;
5677 
5678  nmg_jv(new_v, edge_fus->vp);
5679  edge_fus->vp = new_v;
5680  }
5681 
5682  /* look at all faces around new_v */
5683  vu = BU_LIST_FIRST(vertexuse, &new_v->vu_hd);
5684  while (BU_LIST_NOT_HEAD(vu, &new_v->vu_hd)) {
5685  struct vertexuse *vu_next;
5686  struct faceuse *fu;
5687  struct loopuse *lu;
5688  struct faceuse *bad_fu=(struct faceuse *)NULL;
5689  int bad_loop=0;
5690 
5691  NMG_CK_VERTEXUSE(vu);
5692 
5693  vu_next = BU_LIST_PNEXT(vertexuse, &vu->l);
5694 
5695  if (*vu->up.magic_p != NMG_EDGEUSE_MAGIC) {
5696  vu = vu_next;
5697  continue;
5698  }
5699 
5700  fu = nmg_find_fu_of_vu(vu);
5701  NMG_CK_FACEUSE(fu);
5702 
5703  /* look at all loops in these faces */
5704  lu = BU_LIST_FIRST(loopuse, &fu->lu_hd);
5705  while (BU_LIST_NOT_HEAD(lu, &fu->lu_hd)) {
5706  struct loopuse *lu_next;
5707  struct edgeuse *eu;
5708 
5709  NMG_CK_LOOPUSE(lu);
5710 
5711  lu_next = BU_LIST_PNEXT(loopuse, &lu->l);
5712 
5713  eu = BU_LIST_FIRST(edgeuse, &lu->down_hd);
5714  while (BU_LIST_NOT_HEAD(eu, &lu->down_hd)) {
5715  struct edgeuse *eu_next;
5716 
5717  NMG_CK_EDGEUSE(eu);
5718 
5719  eu_next = BU_LIST_PNEXT(edgeuse, &eu->l);
5720 
5721  /* kill edges that are to/from same vertex */
5722  if (eu->vu_p->v_p == eu->eumate_p->vu_p->v_p) {
5723  while ((vu_next == eu->vu_p || vu_next == eu->eumate_p->vu_p) &&
5724  BU_LIST_NOT_HEAD(vu_next, &new_v->vu_hd))
5725  vu_next = BU_LIST_PNEXT(vertexuse, &vu_next->l);
5726  while ((eu_next == eu || eu_next == eu->eumate_p) &&
5727  BU_LIST_NOT_HEAD(eu_next, &lu->down_hd))
5728  eu_next = BU_LIST_PNEXT(edgeuse, &eu_next->l);
5729 
5730  if (RTG.NMG_debug & DEBUG_BASIC)
5731  bu_log("\tkilling eu %p (%p)\n", (void *)eu, (void *)eu->eumate_p);
5732 
5733  bad_loop = nmg_keu(eu);
5734  }
5735  /* kill edges with length less than tol->dist */
5736  else if (bn_pt3_pt3_equal(eu->vu_p->v_p->vg_p->coord, eu->eumate_p->vu_p->v_p->vg_p->coord, tol)) {
5737  struct edgeuse *prev_eu;
5738 
5739  prev_eu = BU_LIST_PPREV_CIRC(edgeuse, &eu->l);
5740  NMG_CK_EDGEUSE(prev_eu);
5741 
5742  prev_eu->eumate_p->vu_p->v_p = eu->eumate_p->vu_p->v_p;
5743 
5744  while ((vu_next == eu->vu_p || vu_next == eu->eumate_p->vu_p) &&
5745  BU_LIST_NOT_HEAD(vu_next, &new_v->vu_hd))
5746  vu_next = BU_LIST_PNEXT(vertexuse, &vu_next->l);
5747  while ((eu_next == eu || eu_next == eu->eumate_p) &&
5748  BU_LIST_NOT_HEAD(eu_next, &lu->down_hd))
5749  eu_next = BU_LIST_PNEXT(edgeuse, &eu_next->l);
5750 
5751  if (RTG.NMG_debug & DEBUG_BASIC)
5752  bu_log("\tkilling eu %p (%p)\n", (void *)eu, (void *)eu->eumate_p);
5753 
5754  bad_loop = nmg_keu(eu);
5755  }
5756 
5757  if (bad_loop) {
5758  /* emptied a loop, so kill it */
5759  while ((lu_next == lu || lu_next == lu->lumate_p) &&
5760  BU_LIST_NOT_HEAD(lu_next, &fu->lu_hd))
5761  lu_next = BU_LIST_PNEXT(loopuse, &lu_next->l);
5762 
5763  bad_fu = nmg_find_fu_of_lu(lu);
5764  if (!nmg_klu(lu))
5765  bad_fu = (struct faceuse *)NULL;
5766 
5767  break;
5768  }
5769 
5770  eu = eu_next;
5771  }
5772  if (bad_fu) {
5773  /* emptied a faceuse, so kill it */
5774  if (nmg_kfu(bad_fu)) {
5775  /* I can't believe I emptied the whole thing!! */
5776  bu_log("nmg_remove_short_eus_inter: nmg_kfu emptied shell!\n");
5777  break;
5778  }
5779  }
5780  lu = lu_next;
5781  }
5782 
5783  vu = vu_next;
5784  }
5785 }
5786 
5787 
5788 /**
5789  * Eliminates adjacent intersect_fus structs with collinear edges
5790  */
5791 HIDDEN void
5792 nmg_simplify_inter(const struct vertex *new_v, struct bu_ptbl *int_faces, const struct bn_tol *tol)
5793 {
5794  int edge_no=0;
5795  int next_edge_no;
5796 
5797  if (RTG.NMG_debug & DEBUG_BASIC)
5798  bu_log("nmg_simplify_inter(new_v=%p (%f %f %f), int_faces=%p)\n",
5799  (void *)new_v, V3ARGS(new_v->vg_p->coord), (void *)int_faces);
5800 
5801  NMG_CK_VERTEX(new_v);
5802  BU_CK_PTBL(int_faces);
5803  BN_CK_TOL(tol);
5804 
5805  while (BU_PTBL_END(int_faces) > 1 && edge_no < BU_PTBL_END(int_faces)) {
5806  struct intersect_fus *i_fus;
5807  struct intersect_fus *j_fus;
5808 
5809  /* get two consecutive structures */
5810  i_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, edge_no);
5811  next_edge_no = edge_no+1;
5812  if (next_edge_no == BU_PTBL_END(int_faces))
5813  next_edge_no = 0;
5814  j_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, next_edge_no);
5815 
5816  /* skip open space */
5817  if ((i_fus->free_edge || j_fus->free_edge) && next_edge_no == 0) {
5818  edge_no++;
5819  continue;
5820  }
5821 
5822  /* Don't fuse free edges together */
5823  if (i_fus->free_edge && j_fus->free_edge) {
5824  edge_no++;
5825  continue;
5826  }
5827 
5828  /* if either vertex or edgeuse is null, skip */
5829  if (i_fus->vp == NULL || j_fus->vp == NULL ||
5830  i_fus->eu == NULL || j_fus->eu == NULL) {
5831  edge_no++;
5832  continue;
5833  }
5834 
5835  /* If either vertex is new_v, skip */
5836  if (i_fus->vp == new_v || j_fus->vp == new_v) {
5837  edge_no++;
5838  continue;
5839  }
5840 
5841  NMG_CK_VERTEX(i_fus->vp);
5842  NMG_CK_VERTEX(j_fus->vp);
5843  NMG_CK_EDGEUSE(i_fus->eu);
5844  NMG_CK_EDGEUSE(j_fus->eu);
5845 
5846  /* if the two vertices are within tolerance,
5847  * fuse them
5848  */
5849  if (i_fus->vp == j_fus->vp) {
5850  nmg_fuse_inters(i_fus, j_fus, int_faces, tol);
5851  continue;
5852  } else if (bn_pt3_pt3_equal(i_fus->vp->vg_p->coord, j_fus->vp->vg_p->coord, tol)) {
5853  nmg_jv(i_fus->vp, j_fus->vp);
5854  nmg_fuse_inters(i_fus, j_fus, int_faces, tol);
5855  continue;
5856  } else if (bn_3pts_collinear(i_fus->vp->vg_p->coord, j_fus->vp->vg_p->coord, new_v->vg_p->coord, tol)) {
5857  fastf_t i_dist, j_dist;
5858  vect_t i_dist_to_new_v, j_dist_to_new_v;
5859 
5860  /* all three points are collinear,
5861  * may need to split edges
5862  */
5863 
5864  VSUB2(i_dist_to_new_v, new_v->vg_p->coord, i_fus->vp->vg_p->coord);
5865  VSUB2(j_dist_to_new_v, new_v->vg_p->coord, j_fus->vp->vg_p->coord);
5866 
5867  if (VDOT(i_dist_to_new_v, j_dist_to_new_v) < 0.0) {
5868  /* points are collinear with new_v, but in opposite directions */
5869  edge_no++;
5870  continue;
5871  }
5872 
5873  i_dist = MAGSQ(i_dist_to_new_v);
5874  j_dist = MAGSQ(j_dist_to_new_v);
5875 
5876  if (i_dist < tol->dist_sq || j_dist < tol->dist_sq)
5877  bu_bomb("nmg_simplify_inter: vertex within tolerance of new_v\n");
5878 
5879  if (RTG.NMG_debug & DEBUG_BASIC)
5880  bu_log("\tCollinear vertices %p, %p, and %p\n",
5881  (void *)new_v, (void *)i_fus->vp, (void *)j_fus->vp);
5882 
5883  if (i_dist > j_dist && j_dist > tol->dist_sq) {
5884  /* j point is closer to new_v than i point
5885  * split edge at j point
5886  */
5887 
5888  if (RTG.NMG_debug & DEBUG_BASIC)
5889  bu_log("\tSplitting i_fus->eu %p at vertex %p\n",
5890  (void *)i_fus->eu, (void *)j_fus->vp);
5891 
5892  (void)nmg_esplit(j_fus->vp, i_fus->eu, 0);
5893  i_fus->vp = j_fus->vp;
5894  nmg_fuse_inters(i_fus, j_fus, int_faces, tol);
5895 
5896  continue;
5897  } else if (j_dist > i_dist && i_dist > tol->dist_sq) {
5898  /* i point is closer to new_v than j point
5899  * split edge at i point
5900  */
5901 
5902  if (RTG.NMG_debug & DEBUG_BASIC)
5903  bu_log("\tSplitting j_fus->eu %p at vertex %p\n",
5904  (void *)j_fus->eu, (void *)i_fus->vp);
5905 
5906  (void)nmg_esplit(i_fus->vp, j_fus->eu, 0);
5907  nmg_fuse_inters(i_fus, j_fus, int_faces, tol);
5908  continue;
5909  }
5910  }
5911  edge_no++;
5912  }
5913  if (RTG.NMG_debug & DEBUG_BASIC) {
5914  bu_log("\nAfter nmg_simplify_inter:\n");
5915  nmg_pr_inter(new_v, int_faces);
5916  }
5917 }
5918 
5919 
5920 /**
5921  * Make new faces around vertex new_v using info in the table of
5922  * intersect_fu structures. Each structure contains a vertex on an
5923  * edge departing new_v. Vertices from two consecutive edges are
5924  * combined with new_v to form triangular faces around new_v
5925  */
5926 void
5927 nmg_make_faces_at_vert(struct vertex *new_v, struct bu_ptbl *int_faces, const struct bn_tol *tol)
5928 {
5929  struct loopuse *old_lu;
5930  int edge_no=0, next_edge_no;
5931 
5932  if (RTG.NMG_debug & DEBUG_BASIC)
5933  bu_log("nmg_make_faces_at_vert(%p, %ld intersect_fus structs)\n",
5934  (void *)new_v, BU_PTBL_END(int_faces));
5935 
5936  NMG_CK_VERTEX(new_v);
5937  BU_CK_PTBL(int_faces);
5938  BN_CK_TOL(tol);
5939 
5940  if (BU_PTBL_END(int_faces) == 1) {
5941  struct intersect_fus *i_fus;
5942 
5943  /* only one intersect point is left, move new_v to it
5944  * and don't make any faces
5945  */
5946  i_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, 0);
5947  if (i_fus->vp) {
5948  i_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, 0);
5949 
5950  VMOVE(new_v->vg_p->coord, i_fus->vp->vg_p->coord);
5951  nmg_jv(new_v, i_fus->vp);
5952  }
5953  return;
5954  }
5955 
5956  if (BU_PTBL_END(int_faces) == 2) {
5957  struct intersect_fus *i_fus, *j_fus;
5958  point_t center_pt;
5959 
5960  /* only two intersect points left, if they are not on free edges,
5961  * move new_v to the center of the connecting line. No new faces needed
5962  */
5963  i_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, 0);
5964  j_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, 1);
5965 
5966  if (i_fus->vp && j_fus->vp && !i_fus->free_edge && !j_fus->free_edge) {
5967  VCOMB2(center_pt, 0.5, i_fus->vp->vg_p->coord, 0.5, j_fus->vp->vg_p->coord);
5968  VMOVE(new_v->vg_p->coord, center_pt);
5969  }
5970  return;
5971  }
5972 
5973  /* Need to make new faces.
5974  * loop around the vertex, looking at
5975  * pairs of adjacent edges and deciding
5976  * if a new face needs to be constructed
5977  * from the two intersect vertices and new_v
5978  */
5979  while (edge_no < BU_PTBL_END(int_faces)) {
5980  struct intersect_fus *i_fus;
5981  struct intersect_fus *j_fus;
5982  struct vertexuse *vu1, *vu2;
5983  struct edgeuse *eu;
5984  struct loopuse *lu;
5985  struct loopuse *new_lu;
5986  struct faceuse *new_fu;
5987  struct faceuse *fu;
5988 
5989  /* get two consecutive structures */
5990  i_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, edge_no);
5991  next_edge_no = edge_no+1;
5992  if (next_edge_no == BU_PTBL_END(int_faces))
5993  next_edge_no = 0;
5994  j_fus = (struct intersect_fus *)BU_PTBL_GET(int_faces, next_edge_no);
5995 
5996  /* Don't construct a new face across open space */
5997  if ((i_fus->free_edge || j_fus->free_edge) && next_edge_no == 0) {
5998  edge_no++;
5999  continue;
6000  }
6001 
6002  /* if the two vertices are the same, no face needed */
6003  if (i_fus->vp == j_fus->vp) {
6004  edge_no++;
6005  continue;
6006  }
6007 
6008  /* if either vertex is null, no face needed */
6009  if (i_fus->vp == NULL || j_fus->vp == NULL || i_fus->eu == NULL || j_fus->eu == NULL) {
6010  edge_no++;
6011  continue;
6012  }
6013 
6014  /* Don't make faces with two vertices the same */
6015  if (i_fus->vp == new_v || j_fus->vp == new_v) {
6016  edge_no++;
6017  continue;
6018  }
6019 
6020  NMG_CK_VERTEX(i_fus->vp);
6021  NMG_CK_VERTEX(j_fus->vp);
6022  NMG_CK_EDGEUSE(i_fus->eu);
6023  NMG_CK_EDGEUSE(j_fus->eu);
6024 
6025  /* don't make degenerate faces */
6026  if (bn_3pts_collinear(i_fus->vp->vg_p->coord, j_fus->vp->vg_p->coord, new_v->vg_p->coord, tol)) {
6027  edge_no++;
6028  continue;
6029  }
6030 
6031  /* O.K., here is where we actually start making faces.
6032  * Find uses of the two vertices in the same loopuse
6033  */
6034  old_lu = j_fus->eu->up.lu_p;
6035  vu1 = (struct vertexuse *)NULL;
6036  vu2 = (struct vertexuse *)NULL;
6037  for (BU_LIST_FOR (eu, edgeuse, &old_lu->down_hd)) {
6038  if (eu->vu_p->v_p == i_fus->vp)
6039  vu1 = eu->vu_p;
6040  else if (eu->vu_p->v_p == j_fus->vp)
6041  vu2 = eu->vu_p;
6042  }
6043 
6044  if (vu1 == NULL || vu2 == NULL) {
6045  bu_log("nmg_make_faces_at_vert: ERROR: Can't find loop containing vertices %p and %p\n",
6046  (void *)i_fus->vp, (void *)j_fus->vp);
6047  bu_log("\t(%f %f %f) and (%f %f %f)\n", V3ARGS(i_fus->vp->vg_p->coord), V3ARGS(j_fus->vp->vg_p->coord));
6048  edge_no++;
6049  continue;
6050  }
6051 
6052  /* make sure the two vertices have a third between,
6053  * otherwise, don't cut the loop
6054  */
6055  eu = vu1->up.eu_p;
6056  if (eu->eumate_p->vu_p == vu2) {
6057  edge_no++;
6058  continue;
6059  }
6060  eu = vu2->up.eu_p;
6061  if (eu->eumate_p->vu_p == vu1) {
6062  edge_no++;
6063  continue;
6064  }
6065 
6066  /* cut the face loop across the two vertices */
6067  new_lu = nmg_cut_loop(vu1, vu2);
6068 
6069  /* Fix orientations.
6070  * We will never be cutting an OT_OPPOSITE loop
6071  * so the will always be OT_SAME
6072  */
6073  new_lu->orientation = OT_SAME;
6074  new_lu->lumate_p->orientation = OT_SAME;
6075  old_lu->orientation = OT_SAME;
6076  old_lu->lumate_p->orientation = OT_SAME;
6077 
6078  /* find which loopuse contains new_v
6079  * this will be the one to become a new face
6080  */
6081  lu = NULL;
6082 
6083  /* first check old_lu */
6084  for (BU_LIST_FOR (eu, edgeuse, &old_lu->down_hd)) {
6085  if (eu->vu_p->v_p == new_v) {
6086  lu = old_lu;
6087  break;
6088  }
6089  }
6090 
6091  /* if not found check new_lu */
6092  if (lu == NULL) {
6093  for (BU_LIST_FOR (eu, edgeuse, &new_lu->down_hd)) {
6094  if (eu->vu_p->v_p == new_v) {
6095  lu = old_lu;
6096  break;
6097  }
6098  }
6099  }
6100 
6101  if (lu == NULL) {
6102  fu = old_lu->up.fu_p;
6103  bu_log("nmg_make_faces_at_vert: can't find loop for new face\n");
6104  bu_log("vu1 = %p (%p) vu2 = %p (%p)\n",
6105  (void *)vu1, (void *)vu1->v_p, (void *)vu2, (void *)vu2->v_p);
6106  bu_log("new_v = %p\n", (void *)new_v);
6107  bu_log("old_lu = %p, new_lu = %p\n", (void *)old_lu, (void *)new_lu);
6108  nmg_pr_fu_briefly(fu, (char *)NULL);
6109  bu_bomb("nmg_make_faces_at_vert: can't find loop for new face\n");
6110  }
6111 
6112  /* make the new face from the new loop */
6113  new_fu =