libnmg is separated in trunk now (svn branch delete)

svn:revision:69558
svn:branch:libnmg
svn:account:starseeker

.gitignore

@@ -1,3 +0,0 @@
-/bin
-/lib
-/build

include/nmg.h

@@ -33,15 +33,16 @@
 #ifndef NMG_H
 #define NMG_H
 
+
 #include "common.h"
 
 /* interface headers */
 #include "bu/list.h"
-#include "bu/log.h"
 #include "bu/magic.h"
 #include "bu/ptbl.h"
 #include "vmath.h"
 
+
 #define DEBUG_PL_ANIM   0x00000001	/**< @brief 1 mged: animated evaluation */
 #define DEBUG_PL_SLOW   0x00000002	/**< @brief 2 mged: add delays to animation */
 #define DEBUG_GRAPHCL   0x00000004	/**< @brief 3 mged: graphic classification */
@@ -513,18 +514,6 @@ struct vertexuse_a_cnurb {
     long index;			/**< @brief struct # in this model */
 };
 
-/**
- * storage allocation support
- * OBSOLETE
- */
-#define NMG_GETSTRUCT(p, str) BU_GET(p, struct str)
-
-/**
- * storage de-allocation support
- * OBSOLETE
- */
-#define NMG_FREESTRUCT(ptr, str) BU_PUT(ptr, struct str)
-
 
 /*
  * Macros to create and destroy storage for the NMG data structures.
@@ -537,52 +526,52 @@ struct vertexuse_a_cnurb {
 	NMG_CK_MODEL(_m); (_p)->index = ((_m)->maxindex)++
 
 
-#define GET_MODEL(p, m)             {NMG_GETSTRUCT(p, model); NMG_INCR_INDEX(p, m);}
-#define GET_REGION(p, m)            {NMG_GETSTRUCT(p, nmgregion); NMG_INCR_INDEX(p, m);}
-#define GET_REGION_A(p, m)          {NMG_GETSTRUCT(p, nmgregion_a); NMG_INCR_INDEX(p, m);}
-#define GET_SHELL(p, m)             {NMG_GETSTRUCT(p, shell); NMG_INCR_INDEX(p, m);}
-#define GET_SHELL_A(p, m)           {NMG_GETSTRUCT(p, shell_a); NMG_INCR_INDEX(p, m);}
-#define GET_FACE(p, m)              {NMG_GETSTRUCT(p, face); NMG_INCR_INDEX(p, m);}
-#define GET_FACE_G_PLANE(p, m)      {NMG_GETSTRUCT(p, face_g_plane); NMG_INCR_INDEX(p, m);}
-#define GET_FACE_G_SNURB(p, m)      {NMG_GETSTRUCT(p, face_g_snurb); NMG_INCR_INDEX(p, m);}
-#define GET_FACEUSE(p, m)           {NMG_GETSTRUCT(p, faceuse); NMG_INCR_INDEX(p, m);}
-#define GET_LOOP(p, m)              {NMG_GETSTRUCT(p, loop); NMG_INCR_INDEX(p, m);}
-#define GET_LOOP_G(p, m)            {NMG_GETSTRUCT(p, loop_g); NMG_INCR_INDEX(p, m);}
-#define GET_LOOPUSE(p, m)           {NMG_GETSTRUCT(p, loopuse); NMG_INCR_INDEX(p, m);}
-#define GET_LOOPUSE_A(p, m)         {NMG_GETSTRUCT(p, loopuse_a); NMG_INCR_INDEX(p, m);}
-#define GET_EDGE(p, m)              {NMG_GETSTRUCT(p, edge); NMG_INCR_INDEX(p, m);}
-#define GET_EDGE_G_LSEG(p, m)       {NMG_GETSTRUCT(p, edge_g_lseg); NMG_INCR_INDEX(p, m);}
-#define GET_EDGE_G_CNURB(p, m)      {NMG_GETSTRUCT(p, edge_g_cnurb); NMG_INCR_INDEX(p, m);}
-#define GET_EDGEUSE(p, m)           {NMG_GETSTRUCT(p, edgeuse); NMG_INCR_INDEX(p, m);}
-#define GET_VERTEX(p, m)            {NMG_GETSTRUCT(p, vertex); NMG_INCR_INDEX(p, m);}
-#define GET_VERTEX_G(p, m)          {NMG_GETSTRUCT(p, vertex_g); NMG_INCR_INDEX(p, m);}
-#define GET_VERTEXUSE(p, m)         {NMG_GETSTRUCT(p, vertexuse); NMG_INCR_INDEX(p, m);}
-#define GET_VERTEXUSE_A_PLANE(p, m) {NMG_GETSTRUCT(p, vertexuse_a_plane); NMG_INCR_INDEX(p, m);}
-#define GET_VERTEXUSE_A_CNURB(p, m) {NMG_GETSTRUCT(p, vertexuse_a_cnurb); NMG_INCR_INDEX(p, m);}
+#define GET_MODEL(p, m)             {BU_GET(p, struct model); NMG_INCR_INDEX(p, m);}
+#define GET_REGION(p, m)            {BU_GET(p, struct nmgregion); NMG_INCR_INDEX(p, m);}
+#define GET_REGION_A(p, m)          {BU_GET(p, struct nmgregion_a); NMG_INCR_INDEX(p, m);}
+#define GET_SHELL(p, m)             {BU_GET(p, struct shell); NMG_INCR_INDEX(p, m);}
+#define GET_SHELL_A(p, m)           {BU_GET(p, struct shell_a); NMG_INCR_INDEX(p, m);}
+#define GET_FACE(p, m)              {BU_GET(p, struct face); NMG_INCR_INDEX(p, m);}
+#define GET_FACE_G_PLANE(p, m)      {BU_GET(p, struct face_g_plane); NMG_INCR_INDEX(p, m);}
+#define GET_FACE_G_SNURB(p, m)      {BU_GET(p, struct face_g_snurb); NMG_INCR_INDEX(p, m);}
+#define GET_FACEUSE(p, m)           {BU_GET(p, struct faceuse); NMG_INCR_INDEX(p, m);}
+#define GET_LOOP(p, m)              {BU_GET(p, struct loop); NMG_INCR_INDEX(p, m);}
+#define GET_LOOP_G(p, m)            {BU_GET(p, struct loop_g); NMG_INCR_INDEX(p, m);}
+#define GET_LOOPUSE(p, m)           {BU_GET(p, struct loopuse); NMG_INCR_INDEX(p, m);}
+#define GET_LOOPUSE_A(p, m)         {BU_GET(p, struct loopuse_a); NMG_INCR_INDEX(p, m);}
+#define GET_EDGE(p, m)              {BU_GET(p, struct edge); NMG_INCR_INDEX(p, m);}
+#define GET_EDGE_G_LSEG(p, m)       {BU_GET(p, struct edge_g_lseg); NMG_INCR_INDEX(p, m);}
+#define GET_EDGE_G_CNURB(p, m)      {BU_GET(p, struct edge_g_cnurb); NMG_INCR_INDEX(p, m);}
+#define GET_EDGEUSE(p, m)           {BU_GET(p, struct edgeuse); NMG_INCR_INDEX(p, m);}
+#define GET_VERTEX(p, m)            {BU_GET(p, struct vertex); NMG_INCR_INDEX(p, m);}
+#define GET_VERTEX_G(p, m)          {BU_GET(p, struct vertex_g); NMG_INCR_INDEX(p, m);}
+#define GET_VERTEXUSE(p, m)         {BU_GET(p, struct vertexuse); NMG_INCR_INDEX(p, m);}
+#define GET_VERTEXUSE_A_PLANE(p, m) {BU_GET(p, struct vertexuse_a_plane); NMG_INCR_INDEX(p, m);}
+#define GET_VERTEXUSE_A_CNURB(p, m) {BU_GET(p, struct vertexuse_a_cnurb); NMG_INCR_INDEX(p, m);}
 
 
-#define FREE_MODEL(p)             NMG_FREESTRUCT(p, model)
-#define FREE_REGION(p)            NMG_FREESTRUCT(p, nmgregion)
-#define FREE_REGION_A(p)          NMG_FREESTRUCT(p, nmgregion_a)
-#define FREE_SHELL(p)             NMG_FREESTRUCT(p, shell)
-#define FREE_SHELL_A(p)           NMG_FREESTRUCT(p, shell_a)
-#define FREE_FACE(p)              NMG_FREESTRUCT(p, face)
-#define FREE_FACE_G_PLANE(p)      NMG_FREESTRUCT(p, face_g_plane)
-#define FREE_FACE_G_SNURB(p)      NMG_FREESTRUCT(p, face_g_snurb)
-#define FREE_FACEUSE(p)           NMG_FREESTRUCT(p, faceuse)
-#define FREE_LOOP(p)              NMG_FREESTRUCT(p, loop)
-#define FREE_LOOP_G(p)            NMG_FREESTRUCT(p, loop_g)
-#define FREE_LOOPUSE(p)           NMG_FREESTRUCT(p, loopuse)
-#define FREE_LOOPUSE_A(p)         NMG_FREESTRUCT(p, loopuse_a)
-#define FREE_EDGE(p)              NMG_FREESTRUCT(p, edge)
-#define FREE_EDGE_G_LSEG(p)       NMG_FREESTRUCT(p, edge_g_lseg)
-#define FREE_EDGE_G_CNURB(p)      NMG_FREESTRUCT(p, edge_g_cnurb)
-#define FREE_EDGEUSE(p)           NMG_FREESTRUCT(p, edgeuse)
-#define FREE_VERTEX(p)            NMG_FREESTRUCT(p, vertex)
-#define FREE_VERTEX_G(p)          NMG_FREESTRUCT(p, vertex_g)
-#define FREE_VERTEXUSE(p)         NMG_FREESTRUCT(p, vertexuse)
-#define FREE_VERTEXUSE_A_PLANE(p) NMG_FREESTRUCT(p, vertexuse_a_plane)
-#define FREE_VERTEXUSE_A_CNURB(p) NMG_FREESTRUCT(p, vertexuse_a_cnurb)
+#define FREE_MODEL(p)             BU_PUT(p, struct model)
+#define FREE_REGION(p)            BU_PUT(p, struct nmgregion)
+#define FREE_REGION_A(p)          BU_PUT(p, struct nmgregion_a)
+#define FREE_SHELL(p)             BU_PUT(p, struct shell)
+#define FREE_SHELL_A(p)           BU_PUT(p, struct shell_a)
+#define FREE_FACE(p)              BU_PUT(p, struct face)
+#define FREE_FACE_G_PLANE(p)      BU_PUT(p, struct face_g_plane)
+#define FREE_FACE_G_SNURB(p)      BU_PUT(p, struct face_g_snurb)
+#define FREE_FACEUSE(p)           BU_PUT(p, struct faceuse)
+#define FREE_LOOP(p)              BU_PUT(p, struct loop)
+#define FREE_LOOP_G(p)            BU_PUT(p, struct loop_g)
+#define FREE_LOOPUSE(p)           BU_PUT(p, struct loopuse)
+#define FREE_LOOPUSE_A(p)         BU_PUT(p, struct loopuse_a)
+#define FREE_EDGE(p)              BU_PUT(p, struct edge)
+#define FREE_EDGE_G_LSEG(p)       BU_PUT(p, struct edge_g_lseg)
+#define FREE_EDGE_G_CNURB(p)      BU_PUT(p, struct edge_g_cnurb)
+#define FREE_EDGEUSE(p)           BU_PUT(p, struct edgeuse)
+#define FREE_VERTEX(p)            BU_PUT(p, struct vertex)
+#define FREE_VERTEX_G(p)          BU_PUT(p, struct vertex_g)
+#define FREE_VERTEXUSE(p)         BU_PUT(p, struct vertexuse)
+#define FREE_VERTEXUSE_A_PLANE(p) BU_PUT(p, struct vertexuse_a_plane)
+#define FREE_VERTEXUSE_A_CNURB(p) BU_PUT(p, struct vertexuse_a_cnurb)
 
 
 /**
@@ -779,6 +768,10 @@ struct nmg_visit_handlers {
     void (*vis_vertex_g)(uint32_t *, void *, int);
 };
 
+
+extern int nmg_debug;
+
+
 #endif /* NMG_H */
 
 /** @} */

src/CMakeLists.txt

@@ -38,6 +38,7 @@ set(level_1_dirs
   libbn
   libbg
   libbrep
+  libnmg
   librt
   libwdb
   )

src/libnmg/CMakeLists.txt

@@ -0,0 +1,43 @@
+# Include directories needed by libnmg users
+set(NMG_INCLUDE_DIRS
+  ${BRLCAD_BINARY_DIR}/include
+  ${BRLCAD_SOURCE_DIR}/include
+  ${BU_INCLUDE_DIRS}
+  ${BN_INCLUDE_DIRS}
+  ${OPENNURBS_INCLUDE_DIR}
+  )
+BRLCAD_LIB_INCLUDE_DIRS(nmg NMG_INCLUDE_DIRS "")
+
+
+set(LIBNMG_SOURCES
+  nmg_bool.c
+  nmg_brep.cpp
+  nmg_ck.c
+  nmg_class.c
+  nmg_copy.c
+  nmg_eval.c
+  nmg_extrude.c
+  nmg_fcut.c
+  nmg_fuse.c
+  nmg_index.c
+  nmg_info.c
+  nmg_inter.c
+  nmg_manif.c
+  nmg_mesh.c
+  nmg_mirror.c
+  nmg_misc.c
+  nmg_mk.c
+  nmg_mod.c
+  nmg_pr.c
+  nmg_pt_fu.c
+  nmg_plot.c
+  nmg_tri.c
+  nmg_visit.c
+  )
+
+
+BRLCAD_ADDLIB(libnmg "${LIBNMG_SOURCES}" "libbu;libbn")
+set_target_properties(libnmg PROPERTIES VERSION 20.0.1 SOVERSION 20)
+
+
+CMAKEFILES(nmg_extracted.c)

src/libnmg/nmg_brep.cpp

@@ -0,0 +1,263 @@
+/*                    N M G _ B R E P . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2008-2016 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file nmg_brep.cpp
+ *
+ * b-rep support for NMG
+ *
+ */
+
+#include "common.h"
+
+#include "nmg.h"
+#include "nmg_func.h"
+
+#include "bn.h"
+#include "bg.h"
+#include "brep.h"
+#include "bu/log.h"
+#include "bu/malloc.h"
+
+/**
+ *  Routine to create a planar NURBS surface from 4 points
+ */
+HIDDEN ON_Surface*
+sideSurface(const ON_3dPoint& SW, const ON_3dPoint& SE, const ON_3dPoint& NE, const ON_3dPoint& NW)
+{
+    ON_NurbsSurface *surf = ON_NurbsSurface::New(3, 0, 2, 2, 2, 2);
+    surf->SetCV(0, 0, SW);
+    surf->SetCV(1, 0, SE);
+    surf->SetCV(1, 1, NE);
+    surf->SetCV(0, 1, NW);
+    surf->SetKnot(0, 0, 0.0);
+    surf->SetKnot(0, 1, 1.0);
+    surf->SetKnot(1, 0, 0.0);
+    surf->SetKnot(1, 1, 1.0);
+    return surf;
+}
+
+HIDDEN int
+nmg_brep_face(ON_Brep **b, const struct faceuse *fu, const struct bn_tol *tol, long *brepi) {
+    const struct face_g_plane *fg = fu->f_p->g.plane_p;
+    struct bu_ptbl vert_table;
+    struct vertex **pt;
+    int ret = 0;
+    int pnt_cnt = 0;
+    int pnt_index = 0;
+    vect_t u_axis, v_axis;
+    point_t obr_center;
+    point_t *points_3d = NULL;
+    point_t *points_obr = NULL;
+    struct loopuse *lu;
+    struct edgeuse *eu;
+
+    /* Find out how many points we have, set up any uninitialized ON_Brep vertex
+     * structures, and prepare a map of NMG index values to the point array indices */
+    nmg_tabulate_face_g_verts(&vert_table, fg);
+
+    for (BU_PTBL_FOR(pt, (struct vertex **), &vert_table)) {
+	if (brepi[(*pt)->vg_p->index] == -INT_MAX) {
+	    ON_BrepVertex& vert = (*b)->NewVertex((*pt)->vg_p->coord, SMALL_FASTF);
+	    brepi[(*pt)->vg_p->index] = vert.m_vertex_index;
+	}
+	pnt_cnt++;
+    }
+
+    /* Prepare the 3D obr input array */
+    points_3d = (point_t *)bu_calloc(pnt_cnt + 1, sizeof(point_t), "nmg points");
+    for (BU_PTBL_FOR(pt, (struct vertex **), &vert_table)) {
+	VSET(points_3d[pnt_index], (*pt)->vg_p->coord[0],(*pt)->vg_p->coord[1],(*pt)->vg_p->coord[2]);
+	pnt_index++;
+    }
+    bu_ptbl_free(&vert_table);
+
+
+    /* Calculate the 3D coplanar oriented bounding rectangle (obr) */
+    ret += bg_3d_coplanar_obr(&obr_center, &u_axis, &v_axis, (const point_t *)points_3d, pnt_cnt);
+    if (ret) {
+	bu_log("Failed to get oriented bounding rectangle for NMG faceuse #%lu\n", fu->index);
+	return -1;
+    }
+    bu_free(points_3d, "done with obr 3d point inputs");
+
+    /* Use the obr to define the 3D corner points of the NURBS surface */
+    points_obr = (point_t *)bu_calloc(3 + 1, sizeof(point_t), "points_3d");
+    VADD3(points_obr[2], obr_center, u_axis, v_axis);
+    VSCALE(u_axis, u_axis, -1);
+    VADD3(points_obr[3], obr_center, u_axis, v_axis);
+    VSCALE(v_axis, v_axis, -1);
+    VADD3(points_obr[0], obr_center, u_axis, v_axis);
+    VSCALE(u_axis, u_axis, -1);
+    VADD3(points_obr[1], obr_center, u_axis, v_axis);
+
+    /* We need to orient our surface correctly according to the NMG - using
+     * the openNURBS FlipFace function later does not seem to work very
+     * well. If an outer loop is found in the NMG with a cw orientation,
+     * factor that in in addition to the fu->f_p->flip flag. */
+    int ccw = 0;
+    vect_t vtmp, uv1, uv2, vnormal;
+    point_t center;
+    VADD2(center, points_obr[0], points_obr[1]);
+    VADD2(center, center, points_obr[2]);
+    VADD2(center, center, points_obr[3]);
+    VSCALE(center, center, 0.25);
+    for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
+	if (lu->orientation == OT_SAME && nmg_loop_is_ccw(lu, fg->N, tol) == -1) ccw = -1;
+    }
+    if (ccw != -1) {
+	VSET(vnormal, fg->N[0], fg->N[1], fg->N[2]);
+    } else {
+	VSET(vnormal, -fg->N[0], -fg->N[1], -fg->N[2]);
+    }
+    if (fu->f_p->flip)
+	VSET(vnormal, -vnormal[0], -vnormal[1], -vnormal[2]);
+    VSUB2(uv1, points_obr[0], center);
+    VSUB2(uv2, points_obr[1], center);
+    VCROSS(vtmp, uv1, uv2);
+    if (VDOT(vtmp, vnormal) < 0) {
+	VMOVE(vtmp, points_obr[0]);
+	VMOVE(points_obr[0], points_obr[1]);
+	VMOVE(points_obr[1], vtmp);
+	VMOVE(vtmp, points_obr[3]);
+	VMOVE(points_obr[3], points_obr[2]);
+	VMOVE(points_obr[2], vtmp);
+    }
+
+    /* Now that we've got our points correctly oriented for
+     * the NURBS surface, proceed to create it. */
+    ON_3dPoint p1 = ON_3dPoint(points_obr[0]);
+    ON_3dPoint p2 = ON_3dPoint(points_obr[1]);
+    ON_3dPoint p3 = ON_3dPoint(points_obr[2]);
+    ON_3dPoint p4 = ON_3dPoint(points_obr[3]);
+
+    (*b)->m_S.Append(sideSurface(p1, p2, p3, p4));
+    ON_Surface *surf = (*(*b)->m_S.Last());
+    int surfindex = (*b)->m_S.Count();
+    ON_BrepFace& face = (*b)->NewFace(surfindex - 1);
+
+    // With the surface and the face defined, make
+    // trimming loops and create faces.  To generate UV
+    // coordinates for each from and to for the
+    // edgecurves, the UV origin is defined to be v1,
+    // v1->v2 is defined as the U domain, and v1->v4 is
+    // defined as the V domain.
+    VSUB2(u_axis, points_obr[2], points_obr[1]);
+    VSUB2(v_axis, points_obr[0], points_obr[1]);
+    fastf_t u_axis_dist = MAGNITUDE(u_axis);
+    fastf_t v_axis_dist = MAGNITUDE(v_axis);
+
+    /* Now that we have the surface and the face, add the loops */
+    for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
+	if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC) continue; // loop is a single vertex
+	// Check if this is an inner or outer loop
+	ON_BrepLoop::TYPE looptype = (lu->orientation == OT_SAME) ? ON_BrepLoop::outer : ON_BrepLoop::inner;
+	ON_BrepLoop& loop = (*b)->NewLoop(looptype, face);
+	for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+	    vect_t ev1, ev2;
+	    struct vertex_g *vg1 = eu->vu_p->v_p->vg_p;
+	    struct vertex_g *vg2 = eu->eumate_p->vu_p->v_p->vg_p;
+	    NMG_CK_VERTEX_G(vg1);
+	    NMG_CK_VERTEX_G(vg2);
+	    VMOVE(ev1, vg1->coord);
+	    VMOVE(ev2, vg2->coord);
+	    // Add edge if not already added
+	    if (brepi[eu->e_p->index] == -INT_MAX) {
+		/* always add edges with the small vertex index as from */
+		int vert1 = (vg1->index <= vg2->index) ? brepi[vg1->index] : brepi[vg2->index];
+		int vert2 = (vg1->index > vg2->index) ? brepi[vg1->index] : brepi[vg2->index];
+		// Create and add 3D curve
+		ON_Curve* c3d = new ON_LineCurve((*b)->m_V[vert1].Point(), (*b)->m_V[vert2].Point());
+		c3d->SetDomain(0.0, 1.0);
+		(*b)->m_C3.Append(c3d);
+		// Create and add 3D edge
+		ON_BrepEdge& e = (*b)->NewEdge((*b)->m_V[vert1], (*b)->m_V[vert2] , (*b)->m_C3.Count() - 1);
+		e.m_tolerance = 0.0;
+		brepi[eu->e_p->index] = e.m_edge_index;
+	    }
+	    // Regardless of whether the edge existed as an object, it needs to be added to the trimming loop
+	    ON_3dPoint vg1pt(vg1->coord);
+	    int orientation = ((vg1pt != (*b)->m_V[(*b)->m_E[(int)brepi[eu->e_p->index]].m_vi[0]].Point())) ? 1 : 0;
+	    // Make a 2d trimming curve, create a trim, and add the trim to the loop
+	    vect_t vect1, vect2, u_component, v_component;
+	    double u0, u1, v0, v1;
+	    ON_2dPoint from_uv, to_uv;
+	    VSUB2(vect1, ev1, points_obr[0]);
+	    VSUB2(vect2, ev2, points_obr[0]);
+	    surf->GetDomain(0, &u0, &u1);
+	    surf->GetDomain(1, &v0, &v1);
+	    VPROJECT(vect1, u_axis, u_component, v_component);
+	    from_uv.y = u0 + MAGNITUDE(u_component)/u_axis_dist*(u1-u0);
+	    from_uv.x = v0 + MAGNITUDE(v_component)/v_axis_dist*(v1-v0);
+	    VPROJECT(vect2, u_axis, u_component, v_component);
+	    to_uv.y = u0 + MAGNITUDE(u_component)/u_axis_dist*(u1-u0);
+	    to_uv.x = v0 + MAGNITUDE(v_component)/v_axis_dist*(v1-v0);
+	    ON_Curve* c2d =  new ON_LineCurve(from_uv, to_uv);
+	    c2d->SetDomain(0.0, 1.0);
+	    int c2i = (*b)->m_C2.Count();
+	    (*b)->m_C2.Append(c2d);
+	    ON_BrepTrim& trim = (*b)->NewTrim((*b)->m_E[(int)brepi[eu->e_p->index]], orientation, loop, c2i);
+	    trim.m_type = ON_BrepTrim::mated;
+	    trim.m_tolerance[0] = 0.0;
+	    trim.m_tolerance[1] = 0.0;
+	}
+    }
+
+    bu_free(points_obr, "Done with obr");
+
+    return 0;
+}
+
+
+extern "C" void
+rt_nmg_brep(ON_Brep **b, const struct model *m, const struct bn_tol *tol)
+{
+    struct nmgregion *r;
+    struct shell *s;
+    struct faceuse *fu;
+
+    // Verify NMG
+    NMG_CK_MODEL(m);
+
+    // Both NMG and brep structures re-use components between faces.  In order to track
+    // when the conversion routine has already handled an NMG element, use an array.
+    long *brepi = static_cast<long*>(bu_malloc(m->maxindex * sizeof(long), "rt_nmg_brep: brepi[]"));
+    for (int i = 0; i < m->maxindex; i++) brepi[i] = -INT_MAX;
+
+    // Iterate over all faces in the NMG
+    for (BU_LIST_FOR(r, nmgregion, &m->r_hd)) {
+	for (BU_LIST_FOR(s, shell, &r->s_hd)) {
+	    for (BU_LIST_FOR(fu, faceuse, &s->fu_hd)) {
+		NMG_CK_FACEUSE(fu);
+		if (fu->orientation != OT_SAME) continue;
+		if (nmg_brep_face(b, fu, tol, brepi)) return;
+	    }
+	    (*b)->SetTrimIsoFlags();
+	}
+    }
+}
+
+
+// Local Variables:
+// tab-width: 8
+// mode: C++
+// c-basic-offset: 4
+// indent-tabs-mode: t
+// c-file-style: "stroustrup"
+// End:
+// ex: shiftwidth=4 tabstop=8

src/libnmg/nmg_copy.c

@@ -0,0 +1,775 @@
+/*                      N M G _ C O P Y . C
+ * BRL-CAD
+ *
+ * Copyright (c) 2011-2016 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @addtogroup nmg */
+/** @{ */
+/** @file primitives/nmg/nmg_copy.c
+ *
+ * Support routine for n-Manifold Geometry:
+ * Deep copy of NMG model structure
+ *
+ */
+/** @} */
+
+#include "common.h"
+
+#include "nmg.h"
+#include "nmg_func.h"
+
+#include "bu/log.h"
+#include "bu/malloc.h"
+
+#include <string.h>
+
+
+static struct nmgregion_a *
+nmg_construct_region_a(const struct nmgregion_a *original, void **structArray)
+{
+    struct nmgregion_a *ret;
+
+    BU_GET(ret, struct nmgregion_a);
+
+    ret->magic = NMG_REGION_A_MAGIC;
+
+    VMOVE(ret->min_pt, original->min_pt);
+    VMOVE(ret->max_pt, original->max_pt);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct nmgregion *
+nmg_construct_region(struct model *parent, const struct nmgregion *original, void **structArray)
+{
+    struct nmgregion *ret;
+
+    BU_GET(ret, struct nmgregion);
+
+    ret->l.magic = NMG_REGION_MAGIC;
+    ret->m_p     = parent;
+    ret->ra_p    = (struct nmgregion_a *)NULL;
+
+    BU_LIST_INIT(&ret->s_hd);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    if (original->ra_p != NULL) {
+	const struct nmgregion_a *originalAttributes = original->ra_p;
+	struct nmgregion_a *newAttributes
+	    = (struct nmgregion_a *)structArray[originalAttributes->index];
+
+	if (newAttributes == NULL)
+	    newAttributes = nmg_construct_region_a(originalAttributes, structArray);
+
+	ret->ra_p = newAttributes;
+    }
+
+    return ret;
+}
+
+
+static struct face_g_plane *
+nmg_construct_face_g_plane(const struct face_g_plane *original, void **structArray)
+{
+    struct face_g_plane *ret;
+
+    BU_GET(ret, struct face_g_plane);
+
+    ret->magic = NMG_FACE_G_PLANE_MAGIC;
+
+    BU_LIST_INIT(&ret->f_hd);
+    HMOVE(ret->N, original->N);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct face_g_snurb *
+nmg_construct_face_g_snurb(const struct face_g_snurb *original, void **structArray)
+{
+    struct face_g_snurb *ret;
+
+    BU_GET(ret, struct face_g_snurb);
+
+    ret->l.magic = NMG_FACE_G_SNURB_MAGIC;
+
+    BU_LIST_INIT(&ret->f_hd);
+
+    ret->order[0] = original->order[0];
+    ret->order[1] = original->order[1];
+
+    ret->u.magic  = NMG_KNOT_VECTOR_MAGIC;
+    ret->u.k_size = original->u.k_size;
+    ret->u.knots = (fastf_t *)bu_malloc(ret->u.k_size * sizeof(fastf_t),
+					"nmg_construct_face_g_snurb(): u.knots");
+    memcpy(ret->u.knots, original->u.knots, ret->u.k_size * sizeof(fastf_t));
+    ret->v.magic  = NMG_KNOT_VECTOR_MAGIC;
+    ret->v.k_size = original->v.k_size;
+    ret->v.knots = (fastf_t *)bu_malloc(ret->v.k_size * sizeof(fastf_t),
+					"nmg_construct_face_g_snurb(): v.knots");
+    memcpy(ret->v.knots, original->v.knots, ret->v.k_size * sizeof(fastf_t));
+
+    ret->s_size[0]  = original->s_size[0];
+    ret->s_size[1]  = original->s_size[1];
+    ret->pt_type    = original->pt_type;
+    ret->ctl_points
+	= (fastf_t *)bu_malloc(original->s_size[0] * original->s_size[1] * RT_NURB_EXTRACT_COORDS(ret->pt_type) * sizeof(fastf_t),
+			       "nmg_construct_face_g_snurb(): ctl_points");
+    memcpy(ret->ctl_points, original->ctl_points, original->s_size[0] * original->s_size[1] * RT_NURB_EXTRACT_COORDS(ret->pt_type) * sizeof(fastf_t));
+
+    ret->dir = original->dir;
+    VMOVE(ret->min_pt, original->min_pt);
+    VMOVE(ret->max_pt, original->max_pt);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct face *
+nmg_construct_face(struct faceuse *parent, const struct face *original, void **structArray)
+{
+    struct face *ret;
+
+    BU_GET(ret, struct face);
+
+    ret->l.magic   = NMG_FACE_MAGIC;
+    ret->fu_p      = parent;
+    ret->g.magic_p = (uint32_t *)NULL;
+    ret->flip      = original->flip;
+
+    VMOVE(ret->min_pt, original->min_pt);
+    VMOVE(ret->max_pt, original->max_pt);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    switch (*original->g.magic_p) {
+	case NMG_FACE_G_PLANE_MAGIC:
+	    ret->g.plane_p = (struct face_g_plane *)structArray[original->g.plane_p->index];
+
+	    if (ret->g.plane_p == NULL)
+		ret->g.plane_p = nmg_construct_face_g_plane(original->g.plane_p, structArray);
+
+	    BU_LIST_INSERT(&ret->g.plane_p->f_hd, &ret->l);
+	    break;
+
+	case NMG_FACE_G_SNURB_MAGIC:
+	    ret->g.snurb_p = (struct face_g_snurb *)structArray[original->g.plane_p->index];
+
+	    if (ret->g.snurb_p == NULL)
+		ret->g.snurb_p = nmg_construct_face_g_snurb(original->g.snurb_p, structArray);
+
+	    BU_LIST_INSERT(&ret->g.snurb_p->f_hd, &ret->l);
+    }
+
+    return ret;
+}
+
+
+static struct vertex_g *
+nmg_construct_vertex_g(const struct vertex_g *original, void **structArray)
+{
+    struct vertex_g *ret;
+
+    BU_GET(ret, struct vertex_g);
+
+    ret->magic = NMG_VERTEX_G_MAGIC;
+
+    VMOVE(ret->coord, original->coord);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct vertex *
+nmg_construct_vertex(const struct vertex *original, void **structArray)
+{
+    struct vertex *ret;
+
+    BU_GET(ret, struct vertex);
+
+    ret->magic = NMG_VERTEX_MAGIC;
+
+    BU_LIST_INIT(&ret->vu_hd);
+
+    ret->vg_p               = (struct vertex_g*)NULL;
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    if (original->vg_p != NULL) {
+	ret->vg_p = (struct vertex_g *)structArray[original->vg_p->index];
+
+	if (ret->vg_p == NULL)
+	    ret->vg_p = nmg_construct_vertex_g(original->vg_p, structArray);
+    }
+
+    return ret;
+}
+
+
+static struct vertexuse_a_plane *
+nmg_construct_vertexuse_a_plane(const struct vertexuse_a_plane *original, void **structArray)
+{
+    struct vertexuse_a_plane *ret;
+
+    BU_GET(ret, struct vertexuse_a_plane);
+
+    ret->magic            = NMG_VERTEXUSE_A_PLANE_MAGIC;
+
+    VMOVE(ret->N, original->N);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct vertexuse_a_cnurb *
+nmg_construct_vertexuse_a_cnurb(const struct vertexuse_a_cnurb *original, void **structArray)
+{
+    struct vertexuse_a_cnurb *ret;
+
+    BU_GET(ret, struct vertexuse_a_cnurb);
+
+    ret->magic = NMG_VERTEXUSE_A_CNURB_MAGIC;
+
+    VMOVE(ret->param, original->param);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct vertexuse *
+nmg_construct_vertexuse(void *parent, const struct vertexuse *original, void **structArray)
+{
+    struct vertexuse *ret;
+
+    BU_GET(ret, struct vertexuse);
+
+    ret->l.magic            = NMG_VERTEXUSE_MAGIC;
+    ret->up.magic_p         = (uint32_t *)parent;
+    ret->v_p                = (struct vertex*)NULL;
+    ret->a.magic_p          = NULL;
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    if (original->v_p != NULL) {
+	ret->v_p = (struct vertex *)structArray[original->v_p->index];
+
+	if (ret->v_p == NULL)
+	    ret->v_p = nmg_construct_vertex(original->v_p, structArray);
+
+	BU_LIST_INSERT(&ret->v_p->vu_hd, &(ret->l));
+    }
+
+    if (original->a.magic_p != NULL) {
+	switch (*original->a.magic_p) {
+	    case NMG_VERTEXUSE_A_PLANE_MAGIC:
+		ret->a.plane_p = (struct vertexuse_a_plane *)structArray[original->a.plane_p->index];
+		if (ret->a.plane_p == NULL)
+		    ret->a.plane_p = nmg_construct_vertexuse_a_plane(original->a.plane_p, structArray);
+		break;
+	    case NMG_VERTEXUSE_A_CNURB_MAGIC:
+		ret->a.cnurb_p = (struct vertexuse_a_cnurb *)structArray[original->a.cnurb_p->index];
+		if (ret->a.cnurb_p == NULL)
+		    ret->a.cnurb_p = nmg_construct_vertexuse_a_cnurb(original->a.cnurb_p, structArray);
+		break;
+	    default:
+		/* FIXME: any more cases? any action to take? */
+		break;
+	}
+    }
+
+    return ret;
+}
+
+
+static struct edge *
+nmg_construct_edge(struct edgeuse *parent, const struct edge *original, void **structArray)
+{
+    struct edge *ret;
+
+    BU_GET(ret, struct edge);
+
+    ret->magic              = NMG_EDGE_MAGIC;
+    ret->eu_p               = parent;
+    ret->is_real            = original->is_real;
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct edge_g_lseg *
+nmg_construct_edge_g_lseg(const struct edge_g_lseg *original, void **structArray)
+{
+    struct edge_g_lseg *ret;
+
+    BU_GET(ret, struct edge_g_lseg);
+
+    ret->l.magic = NMG_EDGE_G_LSEG_MAGIC;
+
+    BU_LIST_INIT(&ret->eu_hd2);
+
+    VMOVE(ret->e_pt, original->e_pt);
+    VMOVE(ret->e_dir, original->e_dir);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct edge_g_cnurb *
+nmg_construct_edge_g_cnurb(const struct edge_g_cnurb *original, void **structArray)
+{
+    struct edge_g_cnurb *ret;
+
+    BU_GET(ret, struct edge_g_cnurb);
+
+    ret->l.magic = NMG_EDGE_G_CNURB_MAGIC;
+
+    BU_LIST_INIT(&ret->eu_hd2);
+
+    ret->order = original->order;
+
+    ret->k.magic  = NMG_KNOT_VECTOR_MAGIC;
+    ret->k.k_size = original->k.k_size;
+    ret->k.knots = (fastf_t *)bu_malloc(ret->k.k_size * sizeof(fastf_t), "nmg_construct_edge_g_cnurb(): k.knots");
+    memcpy(ret->k.knots, original->k.knots, ret->k.k_size * sizeof(fastf_t));
+
+    ret->c_size     = original->c_size;
+    ret->pt_type    = original->pt_type;
+    ret->ctl_points = (fastf_t *)bu_malloc(ret->c_size * RT_NURB_EXTRACT_COORDS(ret->pt_type) * sizeof(fastf_t),
+					   "nmg_construct_edge_g_cnurb(): ctl_points");
+    memcpy(ret->ctl_points, original->ctl_points, ret->c_size * RT_NURB_EXTRACT_COORDS(ret->pt_type) * sizeof(fastf_t));
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct edgeuse *
+nmg_construct_edgeuse(void *parent, const struct edgeuse *original, void **structArray)
+{
+    struct edgeuse *ret;
+
+    BU_GET(ret, struct edgeuse);
+
+    ret->l.magic = NMG_EDGEUSE_MAGIC;
+
+    BU_LIST_INIT(&ret->l2);
+    ret->l2.magic = NMG_EDGEUSE2_MAGIC;
+
+    ret->up.magic_p         = (uint32_t *)parent;
+    ret->eumate_p           = (struct edgeuse*)NULL;
+    ret->radial_p           = (struct edgeuse*)NULL;
+    ret->e_p                = (struct edge*)NULL;
+    ret->orientation        = original->orientation;
+    ret->vu_p               = (struct vertexuse*) NULL;
+    ret->g.magic_p          = NULL;
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    if (original->eumate_p != NULL) {
+	ret->eumate_p = (struct edgeuse *)structArray[original->eumate_p->index];
+
+	/* because it's tricky to choose the right parent for the mate
+	 * wait until it's created and set eumate_p afterwards
+	 */
+	if (ret->eumate_p != NULL)
+	    ret->eumate_p->eumate_p = ret;
+    }
+
+    if (original->radial_p != NULL) {
+	ret->radial_p = (struct edgeuse *)structArray[original->radial_p->index];
+
+	/* because it's tricky to choose the right parent wait until
+	 * it's created and set it afterwards
+	 */
+	if (ret->radial_p != NULL)
+	    ret->radial_p->radial_p = ret;
+    }
+
+    if (original->e_p != NULL) {
+	ret->e_p = (struct edge *)structArray[original->e_p->index];
+
+	if (ret->e_p == 0)
+	    ret->e_p = nmg_construct_edge(ret, original->e_p, structArray);
+    }
+
+    if (original->vu_p != NULL) {
+	ret->vu_p = (struct vertexuse *)structArray[original->vu_p->index];
+
+	if (ret->vu_p == 0)
+	    ret->vu_p = nmg_construct_vertexuse(ret, original->vu_p, structArray);
+    }
+
+    if (original->g.magic_p != NULL) {
+	switch (*original->g.magic_p) {
+	    case NMG_EDGE_G_LSEG_MAGIC:
+		ret->g.lseg_p = (struct edge_g_lseg *)structArray[original->g.lseg_p->index];
+		if (ret->g.lseg_p == NULL)
+		    ret->g.lseg_p = nmg_construct_edge_g_lseg(original->g.lseg_p, structArray);
+		BU_LIST_INSERT(&ret->g.lseg_p->eu_hd2, &(ret->l2));
+		break;
+	    case NMG_EDGE_G_CNURB_MAGIC:
+		ret->g.cnurb_p = (struct edge_g_cnurb *)structArray[original->g.cnurb_p->index];
+		if (ret->g.cnurb_p == NULL)
+		    ret->g.cnurb_p = nmg_construct_edge_g_cnurb(original->g.cnurb_p, structArray);
+		BU_LIST_INSERT(&ret->g.cnurb_p->eu_hd2, &(ret->l2));
+		break;
+	    default:
+		/* FIXME: any more cases? any action to take? */
+		break;
+	}
+    }
+
+    return ret;
+}
+
+
+static struct loop_g *
+nmg_construct_loop_g(const struct loop_g *original, void **structArray)
+{
+    struct loop_g *ret;
+
+    BU_GET(ret, struct loop_g);
+
+    ret->magic = NMG_LOOP_G_MAGIC;
+
+    VMOVE(ret->min_pt, original->min_pt);
+    VMOVE(ret->max_pt, original->max_pt);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct loop *
+nmg_construct_loop(struct loopuse *parent, const struct loop *original, void **structArray)
+{
+    struct loop *ret;
+
+    BU_GET(ret, struct loop);
+
+    ret->magic              = NMG_LOOP_MAGIC;
+    ret->lu_p               = parent;
+    ret->lg_p               = (struct loop_g*)NULL;
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    if (original->lg_p != NULL) {
+	ret->lg_p = (struct loop_g *)structArray[original->lg_p->index];
+
+	if (ret->lg_p == NULL)
+	    ret->lg_p = nmg_construct_loop_g(original->lg_p, structArray);
+    }
+
+    return ret;
+}
+
+
+static struct loopuse *
+nmg_construct_loopuse(void *parent, const struct loopuse *original, void **structArray)
+{
+    struct loopuse *ret;
+
+    BU_GET(ret, struct loopuse);
+
+    ret->l.magic     = NMG_LOOPUSE_MAGIC;
+    ret->up.magic_p  = (uint32_t *)parent;
+    ret->lumate_p    = (struct loopuse *)NULL;
+    ret->orientation = original->orientation;
+    ret->l_p         = (struct loop *)NULL;
+
+    BU_LIST_INIT(&ret->down_hd);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    if (original->lumate_p != NULL) {
+	ret->lumate_p = (struct loopuse *)structArray[original->lumate_p->index];
+
+	/* because it's tricky to choose the right parent for the mate
+	 * wait until it's created and set eumate_p afterwards
+	 */
+	if (ret->lumate_p != NULL)
+	    ret->lumate_p->lumate_p = ret;
+    }
+
+    if (original->l_p != NULL) {
+	ret->l_p = (struct loop *)structArray[original->l_p->index];
+
+	if (ret->l_p == 0)
+	    ret->l_p = nmg_construct_loop(ret, original->l_p, structArray);
+    }
+
+    switch (BU_LIST_FIRST_MAGIC(&original->down_hd)) {
+	case NMG_VERTEXUSE_MAGIC: {
+	    const struct vertexuse *originalVertexUse = BU_LIST_FIRST(vertexuse, &original->down_hd);
+	    struct vertexuse       *newVertexUse      = (struct vertexuse *)structArray[originalVertexUse->index];
+
+	    if (newVertexUse == NULL)
+		newVertexUse = nmg_construct_vertexuse(ret, originalVertexUse, structArray);
+
+	    BU_LIST_INSERT(&ret->down_hd, &newVertexUse->l);
+	}
+	    break;
+	case NMG_EDGEUSE_MAGIC: {
+	    const struct edgeuse *originalEdgeUse;
+	    for (BU_LIST_FOR(originalEdgeUse, edgeuse, &original->down_hd)) {
+		struct edgeuse *newEdgeUse = (struct edgeuse *)structArray[originalEdgeUse->index];
+		if (newEdgeUse == NULL)
+		    newEdgeUse = nmg_construct_edgeuse(ret, originalEdgeUse, structArray);
+		BU_LIST_INSERT(&ret->down_hd, &newEdgeUse->l);
+	    }
+	}
+	    break;
+	default:
+	    /* FIXME: any more cases? any action to take? */
+	    break;
+    }
+
+    return ret;
+}
+
+
+static struct faceuse *
+nmg_construct_faceuse(struct shell *parent, const struct faceuse *original, void **structArray)
+{
+    struct faceuse       *ret;
+    const struct loopuse *originalLoopUse;
+
+    BU_GET(ret, struct faceuse);
+
+    ret->l.magic     = NMG_FACEUSE_MAGIC;
+    ret->s_p         = parent;
+    ret->fumate_p    = (struct faceuse *)NULL;
+    ret->orientation = original->orientation;
+    ret->outside     = original->outside;
+    ret->f_p         = (struct face *)NULL;
+
+    BU_LIST_INIT(&ret->lu_hd);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    if (original->fumate_p != NULL) {
+	ret->fumate_p = (struct faceuse *)structArray[original->fumate_p->index];
+
+	if (ret->fumate_p == NULL)
+	    ret->fumate_p = nmg_construct_faceuse(parent, original->fumate_p, structArray);
+    }
+
+    if (original->f_p != NULL) {
+	ret->f_p = (struct face *)structArray[original->f_p->index];
+
+	if (ret->f_p == 0)
+	    ret->f_p = nmg_construct_face(ret, original->f_p, structArray);
+    }
+
+    for (BU_LIST_FOR(originalLoopUse, loopuse, &original->lu_hd)) {
+	struct loopuse *newLoopUse = (struct loopuse *)structArray[originalLoopUse->index];
+
+	if (newLoopUse == NULL)
+	    newLoopUse = nmg_construct_loopuse(ret, originalLoopUse, structArray);
+
+	BU_LIST_INSERT(&ret->lu_hd, &newLoopUse->l);
+    }
+
+    return ret;
+}
+
+
+static struct shell_a *
+nmg_construct_shell_a(const struct shell_a *original, void **structArray)
+{
+    struct shell_a *ret;
+
+    BU_GET(ret, struct shell_a);
+
+    ret->magic = NMG_SHELL_A_MAGIC;
+
+    VMOVE(ret->min_pt, original->min_pt);
+    VMOVE(ret->max_pt, original->max_pt);
+
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    return ret;
+}
+
+
+static struct shell *
+nmg_construct_shell(struct nmgregion *parent, const struct shell *original, void **structArray)
+{
+    struct shell         *ret;
+    const struct faceuse *originalFaceUse;
+    const struct loopuse *originalLoopUse;
+    const struct edgeuse *originalEdgeUse;
+
+    BU_GET(ret, struct shell);
+
+    ret->l.magic = NMG_SHELL_MAGIC;
+    ret->r_p     = parent;
+    ret->sa_p    = (struct shell_a *)NULL;
+
+    BU_LIST_INIT(&ret->fu_hd);
+    BU_LIST_INIT(&ret->lu_hd);
+    BU_LIST_INIT(&ret->eu_hd);
+
+    ret->vu_p               = (struct vertexuse *) NULL;
+    ret->index              = original->index;
+    structArray[ret->index] = ret;
+
+    if (original->sa_p != NULL) {
+	const struct shell_a *originalAttributes = original->sa_p;
+	struct shell_a       *newAttributes      = (struct shell_a *)structArray[originalAttributes->index];
+
+	if (newAttributes == NULL)
+	    newAttributes = nmg_construct_shell_a(originalAttributes, structArray);
+
+	ret->sa_p = newAttributes;
+    }
+
+    for (BU_LIST_FOR(originalFaceUse, faceuse, &original->fu_hd)) {
+	struct faceuse *newFaceUse = (struct faceuse *)structArray[originalFaceUse->index];
+
+	if (newFaceUse == NULL)
+	    newFaceUse = nmg_construct_faceuse(ret, originalFaceUse, structArray);
+
+	BU_LIST_INSERT(&ret->fu_hd, &newFaceUse->l);
+    }
+
+    for (BU_LIST_FOR(originalLoopUse, loopuse, &original->lu_hd)) {
+	struct loopuse *newLoopUse = (struct loopuse *)structArray[originalLoopUse->index];
+
+	if (newLoopUse == NULL)
+	    newLoopUse = nmg_construct_loopuse(ret, originalLoopUse, structArray);
+
+	BU_LIST_INSERT(&ret->lu_hd, &newLoopUse->l);
+    }
+
+    for (BU_LIST_FOR(originalEdgeUse, edgeuse, &original->eu_hd)) {
+	struct edgeuse *newEdgeUse = (struct edgeuse *)structArray[originalEdgeUse->index];
+
+	if (newEdgeUse == NULL)
+	    newEdgeUse = nmg_construct_edgeuse(ret, originalEdgeUse, structArray);
+
+	BU_LIST_INSERT(&ret->eu_hd, &newEdgeUse->l);
+    }
+
+    if (original->vu_p != 0) {
+	ret->vu_p = (struct vertexuse *)structArray[original->vu_p->index];
+
+	if (ret->vu_p == NULL)
+	    ret->vu_p = nmg_construct_vertexuse(ret, original->vu_p, structArray);
+    }
+
+    return ret;
+}
+
+
+/**
+ * Makes a deep copy of a NMG model structure.
+ */
+struct model *
+nmg_clone_model(const struct model *original)
+{
+    struct model           *ret;
+    void *              *structArray;
+    const struct nmgregion *originalRegion;
+    struct bn_tol           tolerance;
+
+    NMG_CK_MODEL(original);
+
+    structArray = (void **)bu_calloc(original->maxindex, sizeof(void *), "nmg_clone_model() structArray");
+
+    ret = nmg_mm();
+    ret->index    = original->index;
+    ret->maxindex = original->maxindex;
+
+    structArray[ret->index] = ret;
+
+    tolerance.magic   = BN_TOL_MAGIC;
+    tolerance.dist    = 0.0005;
+    tolerance.dist_sq = tolerance.dist * tolerance.dist;
+    tolerance.perp    = 1e-6;
+    tolerance.para    = 1 - tolerance.perp;
+
+    for (BU_LIST_FOR(originalRegion, nmgregion, &original->r_hd)) {
+	struct nmgregion *newRegion = (struct nmgregion *)structArray[originalRegion->index];
+
+	if (newRegion == NULL) {
+	    const struct shell *originalShell;
+
+	    newRegion = nmg_construct_region(ret, originalRegion, structArray);
+
+	    for (BU_LIST_FOR(originalShell, shell, &originalRegion->s_hd)) {
+		struct shell *newShell = (struct shell *)structArray[originalShell->index];
+
+		if (newShell == NULL)
+		    newShell = nmg_construct_shell(newRegion, originalShell, structArray);
+
+		BU_LIST_INSERT(&newRegion->s_hd, &newShell->l);
+	    }
+
+	    BU_LIST_INSERT(&ret->r_hd, &newRegion->l);
+	}
+    }
+
+    bu_free(structArray, "nmg_clone_model() structArray");
+
+    return ret;
+}
+
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

src/libnmg/nmg_eval.c

@@ -0,0 +1,657 @@
+/*                      N M G _ E V A L . C
+ * BRL-CAD
+ *
+ * Copyright (c) 1990-2016 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @addtogroup nmg */
+/** @{ */
+/** @file primitives/nmg/nmg_eval.c
+ *
+ * Evaluate boolean operations on NMG objects.
+ *
+ *
+ */
+/** @} */
+
+#include "common.h"
+
+#include "nmg.h"
+#include "nmg_func.h"
+
+#include "bio.h"
+#include "bu/log.h"
+
+#include <stdlib.h>
+#include <math.h>
+
+
+struct nmg_bool_state {
+    struct shell *bs_dest;
+    struct shell *bs_src;
+    int bs_isA;		/* true if A, else doing B */
+    char **bs_classtab;
+    const int *bs_actions;
+    const struct bn_tol *bs_tol;
+};
+
+
+HIDDEN void nmg_eval_shell(struct shell *s, struct nmg_bool_state *bs);
+HIDDEN void nmg_eval_plot(struct nmg_bool_state *bs, int num);
+
+
+#define BACTION_KILL 1
+#define BACTION_RETAIN 2
+#define BACTION_RETAIN_AND_FLIP 3
+
+static const char *nmg_baction_names[] = {
+    "*undefined 0*",
+    "BACTION_KILL",
+    "BACTION_RETAIN",
+    "BACTION_RETAIN_AND_FLIP",
+    "*undefined 4*"
+};
+
+
+#define NMG_CLASS_BAD 8
+static const char *nmg_class_names[] = {
+    "onAinB",
+    "onAonBshared",
+    "onAonBanti",
+    "onAoutB",
+    "inAonB",
+    "onAonBshared",
+    "onAonBanti",
+    "outAonB",
+    "*BAD*CLASS*"
+};
+
+
+/**
+ * Make sure that the lu and fu orientation flags are consistent with
+ * the geometric arrangement of the vertices and the faceuse normal.
+ */
+void
+nmg_ck_lu_orientation(struct loopuse *lu, const struct bn_tol *tolp)
+{
+    struct faceuse *fu;
+    plane_t fu_peqn;
+    plane_t lu_peqn;
+    fastf_t dot;
+
+    NMG_CK_LOOPUSE(lu);
+    fu = lu->up.fu_p;		/* parent had better be faceuse */
+    NMG_CK_FACEUSE(fu);
+
+    NMG_GET_FU_PLANE(fu_peqn, fu);
+    nmg_loop_plane_newell(lu, lu_peqn);
+
+    dot = VDOT(fu_peqn, lu_peqn);
+
+    if (NEAR_ZERO(dot, tolp->perp))
+	return;		/* can't determine geometric orientation */
+
+
+    if (dot < 0.0) {
+	bu_log("nmg_ck_lu_orientation() lu=%p, dot=%g, fu_orient=%s, lu_orient=%s\n",
+	       (void *)lu, dot,
+	       nmg_orientation(fu->orientation),
+	       nmg_orientation(lu->orientation)
+	    );
+	bu_bomb("nmg_ck_lu_orientation() loop orientation flags do not match geometry\n");
+    }
+}
+
+
+/**
+ * Convert an NMG_CLASS_xxx token into a string name.
+ */
+const char *
+nmg_class_name(int nmg_class)
+{
+    if (nmg_class == NMG_CLASS_Unknown) return "Unknown";
+    if (nmg_class < 0 || nmg_class > NMG_CLASS_BAD) nmg_class = NMG_CLASS_BAD;
+    return nmg_class_names[nmg_class];
+}
+
+
+/*
+ * Action Table for Boolean Operations.
+ *
+ * Each table lists what actions are to be taken for topological elements
+ * which have have each kind of classification.
+ *
+ * Actions are listed in this order:
+ * (Aon) onAinB, onAonBshared, onAonBanti-shared, onAoutB,
+ * (Bon) inAonB, onAonBshared, onAonBanti-shared, outAonB
+ */
+static const int subtraction_actions[8] = {
+    BACTION_KILL,
+    BACTION_KILL,		/* shared */
+    BACTION_RETAIN,		/* anti-shared */
+    BACTION_RETAIN,
+
+    BACTION_RETAIN,		/* (formerly BACTION_RETAIN_AND_FLIP) */
+    BACTION_KILL,
+    BACTION_KILL,
+    BACTION_KILL
+};
+
+
+static const int union_actions[8] = {
+    BACTION_KILL,
+    BACTION_RETAIN,		/* shared */
+    BACTION_KILL,		/* anti-shared */
+    BACTION_RETAIN,
+
+    BACTION_KILL,
+    BACTION_KILL,
+    BACTION_KILL,
+    BACTION_RETAIN
+};
+
+
+static const int intersect_actions[8] = {
+    BACTION_RETAIN,
+    BACTION_RETAIN,		/* shared */
+    BACTION_KILL,		/* anti-shared ==> non-manifold result */
+    BACTION_KILL,
+
+    BACTION_RETAIN,
+    BACTION_KILL,
+    BACTION_KILL,
+    BACTION_KILL
+};
+
+
+/**
+ * Evaluate a boolean operation on the two shells "A" and "B", of the
+ * form "answer = A op B".  As input, each element (loop-in-face, wire
+ * loop, wire edge, vertex) in both A and B has been classified as
+ * being "in", "on", or "out" of the other shell.  Using these
+ * classifications, operate on the input shells.  At the end, shell A
+ * contains the resultant object, and shell B is destroyed.
+ *
+ */
+void
+nmg_evaluate_boolean(struct shell *sA, struct shell *sB, int op, char **classlist, const struct bn_tol *tol)
+{
+    int const *actions;
+    struct nmg_bool_state bool_state;
+
+    NMG_CK_SHELL(sA);
+    NMG_CK_SHELL(sB);
+    BN_CK_TOL(tol);
+
+    if (nmg_debug & DEBUG_BOOLEVAL) {
+	bu_log("nmg_evaluate_boolean(sA=%p, sB=%p, op=%d) START\n",
+	       (void *)sA, (void *)sB, op);
+    }
+
+    switch (op) {
+	case NMG_BOOL_SUB:
+	    actions = subtraction_actions;
+	    nmg_invert_shell(sB);	/* FLIP all faceuse normals */
+	    break;
+	case NMG_BOOL_ADD:
+	    actions = union_actions;
+	    break;
+	case NMG_BOOL_ISECT:
+	    actions = intersect_actions;
+	    break;
+	default:
+	    actions = union_actions;	/* shut up lint */
+	    bu_log("ERROR nmg_evaluate_boolean() op=%d.\n", op);
+	    bu_bomb("bad boolean\n");
+    }
+
+    bool_state.bs_dest = sA;
+    bool_state.bs_src = sB;
+    bool_state.bs_classtab = classlist;
+    bool_state.bs_actions = actions;
+    bool_state.bs_tol = tol;
+
+    bool_state.bs_isA = 1;
+    nmg_eval_shell(sA, &bool_state);
+
+    bool_state.bs_isA = 0;
+    nmg_eval_shell(sB, &bool_state);
+
+    if (nmg_debug & DEBUG_BOOLEVAL) {
+	bu_log("nmg_evaluate_boolean(sA=%p, sB=%p, op=%d), evaluations done\n",
+	       (void *)sA, (void *)sB, op);
+    }
+    /* Write sA and sB into separate files, if wanted? */
+
+    /* Move everything left in sB into sA.  sB is killed. */
+    nmg_js(sA, sB, tol);
+
+    /* Plot the result */
+    if (nmg_debug & DEBUG_BOOLEVAL && nmg_debug & DEBUG_PLOTEM) {
+	FILE *fp;
+
+	if ((fp=fopen("bool_ans.plot3", "wb")) == (FILE *)NULL) {
+	    (void)perror("bool_ans.plot3");
+	    bu_bomb("unable to open bool_ans.plot3 for writing");
+	}
+	bu_log("plotting bool_ans.plot3\n");
+	nmg_pl_s(fp, sA);
+	(void)fclose(fp);
+    }
+
+    /* Remove loops/edges/vertices that appear more than once in result */
+    nmg_rm_redundancies(sA, tol);
+}
+
+
+static int nmg_eval_count = 0;	/* debug -- plot file numbering */
+
+
+/**
+ * Given a pointer to some NMG data structure, search the 4
+ * classification lists to determine its classification.  (XXX In the
+ * future, this should be done with one big array).  Then, return the
+ * action code for an item of that classification.
+ */
+HIDDEN int
+nmg_eval_action(uint32_t *ptr, register struct nmg_bool_state *bs)
+{
+    register int ret;
+    register int nmg_class;
+    int index;
+
+    BN_CK_TOL(bs->bs_tol);
+
+    index = nmg_index_of_struct(ptr);
+    if (bs->bs_isA) {
+	if (NMG_INDEX_VALUE(bs->bs_classtab[NMG_CLASS_AinB], index)) {
+	    nmg_class = NMG_CLASS_AinB;
+	    ret = bs->bs_actions[NMG_CLASS_AinB];
+	    goto out;
+	}
+	if (NMG_INDEX_VALUE(bs->bs_classtab[NMG_CLASS_AonBshared], index)) {
+	    nmg_class = NMG_CLASS_AonBshared;
+	    ret = bs->bs_actions[NMG_CLASS_AonBshared];
+	    goto out;
+	}
+	if (NMG_INDEX_VALUE(bs->bs_classtab[NMG_CLASS_AonBanti], index)) {
+	    nmg_class = NMG_CLASS_AonBanti;
+	    ret = bs->bs_actions[NMG_CLASS_AonBanti];
+	    goto out;
+	}
+	if (NMG_INDEX_VALUE(bs->bs_classtab[NMG_CLASS_AoutB], index)) {
+	    nmg_class = NMG_CLASS_AoutB;
+	    ret = bs->bs_actions[NMG_CLASS_AoutB];
+	    goto out;
+	}
+	bu_log("nmg_eval_action(ptr=%p) %s has no A classification, retaining\n",
+	       (void *)ptr, bu_identify_magic(*((uint32_t *)ptr)));
+	nmg_class = NMG_CLASS_BAD;
+	ret = BACTION_RETAIN;
+	goto out;
+    }
+
+    /* is B */
+    if (NMG_INDEX_VALUE(bs->bs_classtab[NMG_CLASS_BinA], index)) {
+	nmg_class = NMG_CLASS_BinA;
+	ret = bs->bs_actions[NMG_CLASS_BinA];
+	goto out;
+    }
+    if (NMG_INDEX_VALUE(bs->bs_classtab[NMG_CLASS_BonAshared], index)) {
+	nmg_class = NMG_CLASS_BonAshared;
+	ret = bs->bs_actions[NMG_CLASS_BonAshared];
+	goto out;
+    }
+    if (NMG_INDEX_VALUE(bs->bs_classtab[NMG_CLASS_BonAanti], index)) {
+	nmg_class = NMG_CLASS_BonAanti;
+	ret = bs->bs_actions[NMG_CLASS_BonAanti];
+	goto out;
+    }
+    if (NMG_INDEX_VALUE(bs->bs_classtab[NMG_CLASS_BoutA], index)) {
+	nmg_class = NMG_CLASS_BoutA;
+	ret = bs->bs_actions[NMG_CLASS_BoutA];
+	goto out;
+    }
+    bu_log("nmg_eval_action(ptr=%p) %s has no B classification, retaining\n",
+	   (void *)ptr, bu_identify_magic(*((uint32_t *)ptr)));
+    nmg_class = NMG_CLASS_BAD;
+    ret = BACTION_RETAIN;
+out:
+    if (nmg_debug & DEBUG_BOOLEVAL) {
+	bu_log("nmg_eval_action(ptr=%p) index=%d %s %s %s %s\n",
+	       (void *)ptr, index,
+	       bs->bs_isA ? "A" : "B",
+	       bu_identify_magic(*((uint32_t *)ptr)),
+	       nmg_class_name(nmg_class),
+	       nmg_baction_names[ret]);
+    }
+    return ret;
+}
+
+
+/**
+ * Make a life-and-death decision on every element of a shell.
+ * Descend the "great chain of being" from the face to loop to edge to
+ * vertex, saving or demoting along the way.
+ *
+ * Note that there is no moving of items from one shell to another.
+ */
+HIDDEN void
+nmg_eval_shell(register struct shell *s, struct nmg_bool_state *bs)
+{
+    struct faceuse *fu;
+    struct faceuse *nextfu;
+    struct loopuse *lu;
+    struct loopuse *nextlu;
+    struct edgeuse *eu;
+    struct edgeuse *nexteu;
+    struct vertexuse *vu;
+    int loops_retained;
+
+    NMG_CK_SHELL(s);
+    BN_CK_TOL(bs->bs_tol);
+
+    if (nmg_debug & DEBUG_VERIFY)
+	nmg_vshell(&s->r_p->s_hd, s->r_p);
+
+    /*
+     * For each face in the shell, process all the loops in the face,
+     * and then handle the face and all loops as a unit.
+     */
+    nmg_eval_plot(bs, nmg_eval_count++);	/* debug */
+    fu = BU_LIST_FIRST(faceuse, &s->fu_hd);
+    while (BU_LIST_NOT_HEAD(fu, &s->fu_hd)) {
+	NMG_CK_FACEUSE(fu);
+	nextfu = BU_LIST_PNEXT(faceuse, fu);
+
+	/* Faceuse mates will be handled at same time as OT_SAME fu */
+	if (fu->orientation != OT_SAME) {
+	    fu = nextfu;
+	    continue;
+	}
+	if (fu->fumate_p == nextfu)
+	    nextfu = BU_LIST_PNEXT(faceuse, nextfu);
+
+	/* Consider this face */
+	NMG_CK_FACE(fu->f_p);
+
+	loops_retained = 0;
+	lu = BU_LIST_FIRST(loopuse, &fu->lu_hd);
+	while (BU_LIST_NOT_HEAD(lu, &fu->lu_hd)) {
+	    NMG_CK_LOOPUSE(lu);
+	    nextlu = BU_LIST_PNEXT(loopuse, lu);
+	    if (lu->lumate_p == nextlu)
+		nextlu = BU_LIST_PNEXT(loopuse, nextlu);
+
+	    NMG_CK_LOOP(lu->l_p);
+	    nmg_ck_lu_orientation(lu, bs->bs_tol);
+	    switch (nmg_eval_action(&lu->l_p->magic, bs)) {
+		case BACTION_KILL:
+		    /* Kill by demoting loop to edges */
+		    if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+			/* loop of single vertex */
+			(void)nmg_klu(lu);
+		    } else if (nmg_demote_lu(lu) == 0) {
+			nmg_eval_plot(bs, nmg_eval_count++);	/* debug */
+		    }
+		    lu = nextlu;
+		    continue;
+		case BACTION_RETAIN:
+		    loops_retained++;
+		    break;
+		default:
+		    bu_bomb("nmg_eval_shell() bad BACTION\n");
+	    }
+	    lu = nextlu;
+	}
+
+	if (nmg_debug & DEBUG_BOOLEVAL)
+	    bu_log("faceuse %p loops retained=%d\n",
+		   (void *)fu, loops_retained);
+	if (nmg_debug & DEBUG_VERIFY)
+	    nmg_vshell(&s->r_p->s_hd, s->r_p);
+
+	/*
+	 * Here, faceuse will have 0 or more loopuses still in it.
+	 * Decide the fate of the face;  if the face dies,
+	 * then any remaining loops, edges, etc., will die too.
+	 */
+	if (BU_LIST_IS_EMPTY(&fu->lu_hd)) {
+	    if (loops_retained) bu_bomb("nmg_eval_shell() empty faceuse with retained loops?\n");
+	    /* faceuse is empty, face & mate die */
+	    if (nmg_debug & DEBUG_BOOLEVAL)
+		bu_log("faceuse %p empty, kill\n", (void *)fu);
+	    nmg_kfu(fu);	/* kill face & mate, dequeue from shell */
+	    if (nmg_debug & DEBUG_VERIFY)
+		nmg_vshell(&s->r_p->s_hd, s->r_p);
+	    nmg_eval_plot(bs, nmg_eval_count++);	/* debug */
+	    fu = nextfu;
+	    continue;
+	}
+
+	if (loops_retained <= 0) {
+	    nmg_pr_fu(fu, (char *)NULL);
+	    bu_bomb("nmg_eval_shell() non-empty faceuse, no loops retained?\n");
+	}
+	fu = nextfu;
+    }
+    if (nmg_debug & DEBUG_VERIFY)
+	nmg_vshell(&s->r_p->s_hd, s->r_p);
+
+    /*
+     * For each loop in the shell, process.
+     * Each loop is either a wire-loop, or a vertex-with-self-loop.
+     * Only consider wire loops here.
+     */
+    nmg_eval_plot(bs, nmg_eval_count++);	/* debug */
+    lu = BU_LIST_FIRST(loopuse, &s->lu_hd);
+    while (BU_LIST_NOT_HEAD(lu, &s->lu_hd)) {
+	NMG_CK_LOOPUSE(lu);
+	nextlu = BU_LIST_PNEXT(loopuse, lu);
+	if (lu->lumate_p == nextlu)
+	    nextlu = BU_LIST_PNEXT(loopuse, nextlu);
+
+	if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+	    /* ignore vertex-with-self-loop */
+	    lu = nextlu;
+	    continue;
+	}
+	NMG_CK_LOOP(lu->l_p);
+	switch (nmg_eval_action(&lu->l_p->magic, bs)) {
+	    case BACTION_KILL:
+		/* Demote the loopuse into wire edges */
+		/* kill loop & mate */
+		if (nmg_demote_lu(lu) == 0)
+		    nmg_eval_plot(bs, nmg_eval_count++);	/* debug */
+		lu = nextlu;
+		continue;
+	    case BACTION_RETAIN:
+		break;
+	    default:
+		bu_bomb("nmg_eval_shell() bad BACTION\n");
+	}
+	lu = nextlu;
+    }
+    if (nmg_debug & DEBUG_VERIFY)
+	nmg_vshell(&s->r_p->s_hd, s->r_p);
+
+    /*
+     * For each wire-edge in the shell, ...
+     */
+    nmg_eval_plot(bs, nmg_eval_count++);	/* debug */
+    eu = BU_LIST_FIRST(edgeuse, &s->eu_hd);
+    while (BU_LIST_NOT_HEAD(eu, &s->eu_hd)) {
+	NMG_CK_EDGEUSE(eu);
+	nexteu = BU_LIST_PNEXT(edgeuse, eu);	/* may be head */
+	if (eu->eumate_p == nexteu)
+	    nexteu = BU_LIST_PNEXT(edgeuse, nexteu);
+
+	/* Consider this edge */
+	NMG_CK_EDGE(eu->e_p);
+	switch (nmg_eval_action(&eu->e_p->magic, bs)) {
+	    case BACTION_KILL:
+		/* Demote the edegeuse (and mate) into vertices */
+		if (nmg_demote_eu(eu) == 0)
+		    nmg_eval_plot(bs, nmg_eval_count++);	/* debug */
+		eu = nexteu;
+		continue;
+	    case BACTION_RETAIN:
+		break;
+	    default:
+		bu_bomb("nmg_eval_shell() bad BACTION\n");
+	}
+	eu = nexteu;
+    }
+
+    /*
+     * For each lone vertex-with-self-loop, process.
+     * Note that these are intermixed in the loop list.
+     * Each loop is either a wire-loop, or a vertex-with-self-loop.
+     * Only consider cases of vertex-with-self-loop here.
+     *
+     * This case has to be handled separately, because a wire-loop
+     * may be demoted to a set of wire-edges above, some of which
+     * may be retained.  The non-retained wire-edges may in turn
+     * be demoted into vertex-with-self-loop objects above,
+     * which will be processed here.
+     */
+    nmg_eval_plot(bs, nmg_eval_count++);	/* debug */
+    lu = BU_LIST_FIRST(loopuse, &s->lu_hd);
+    while (BU_LIST_NOT_HEAD(lu, &s->lu_hd)) {
+	NMG_CK_LOOPUSE(lu);
+	nextlu = BU_LIST_PNEXT(loopuse, lu);
+
+	if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_VERTEXUSE_MAGIC) {
+	    /* ignore any remaining wire-loops */
+	    lu = nextlu;
+	    continue;
+	}
+	if (nextlu == lu->lumate_p)
+	    nextlu = BU_LIST_PNEXT(loopuse, nextlu);
+	vu = BU_LIST_PNEXT(vertexuse, &lu->down_hd);
+	NMG_CK_VERTEXUSE(vu);
+	NMG_CK_VERTEX(vu->v_p);
+	switch (nmg_eval_action(&vu->v_p->magic, bs)) {
+	    case BACTION_KILL:
+		/* Eliminate the loopuse, and mate */
+		nmg_klu(lu);
+		lu = nextlu;
+		continue;
+	    case BACTION_RETAIN:
+		break;
+	    default:
+		bu_bomb("nmg_eval_shell() bad BACTION\n");
+	}
+	lu = nextlu;
+    }
+    if (nmg_debug & DEBUG_VERIFY)
+	nmg_vshell(&s->r_p->s_hd, s->r_p);
+
+    /*
+     * Final case:  shell of a single vertexuse
+     */
+    vu = s->vu_p;
+    if (vu) {
+	NMG_CK_VERTEXUSE(vu);
+	NMG_CK_VERTEX(vu->v_p);
+	switch (nmg_eval_action(&vu->v_p->magic, bs)) {
+	    case BACTION_KILL:
+		nmg_kvu(vu);
+		nmg_eval_plot(bs, nmg_eval_count++);	/* debug */
+		s->vu_p = (struct vertexuse *)0;	/* sanity */
+		break;
+	    case BACTION_RETAIN:
+		break;
+	    default:
+		bu_bomb("nmg_eval_shell() bad BACTION\n");
+	}
+    }
+    if (nmg_debug & DEBUG_VERIFY)
+	nmg_vshell(&s->r_p->s_hd, s->r_p);
+    nmg_eval_plot(bs, nmg_eval_count++);	/* debug */
+}
+
+
+/**
+ * Called from nmg_eval_shell
+ *
+ * Located here because definition of nmg_bool_state is local to this
+ * module.
+ */
+HIDDEN void
+nmg_eval_plot(struct nmg_bool_state *bs, int num)
+{
+    FILE *fp;
+    char fname[128];
+    int do_plot = 0;
+    int do_anim = 0;
+
+    if (nmg_debug & DEBUG_BOOLEVAL && nmg_debug & DEBUG_PLOTEM)
+	do_plot = 1;
+    if (nmg_debug & DEBUG_PL_ANIM) do_anim = 1;
+
+    if (!do_plot && !do_anim) return;
+
+    BN_CK_TOL(bs->bs_tol);
+
+    if (do_plot) {
+	sprintf(fname, "nmg_eval%d.plot3", num);
+	if ((fp = fopen(fname, "wb")) == NULL) {
+	    perror(fname);
+	    return;
+	}
+	bu_log("Plotting %s\n", fname);
+
+	nmg_pl_s(fp, bs->bs_dest);
+	nmg_pl_s(fp, bs->bs_src);
+
+	fclose(fp);
+    }
+
+    if (do_anim) {
+	struct bn_vlblock *vbp;
+
+	vbp = rt_vlblock_init();
+
+	nmg_vlblock_s(vbp, bs->bs_dest, 0);
+	nmg_vlblock_s(vbp, bs->bs_src, 0);
+
+	/* Cause animation of boolean operation as it proceeds! */
+	if (nmg_vlblock_anim_upcall) {
+	    /* if requested, delay 1/4 second */
+	    /* need to cast nmg_vlblock_anim_upcall pointer for actual use as a function */
+	    void (*cfp)(struct bn_vlblock *, int, int);
+	    cfp = (void (*)(struct bn_vlblock *, int, int))nmg_vlblock_anim_upcall;
+	    cfp(vbp,
+		(nmg_debug&DEBUG_PL_SLOW) ? 250000 : 0,
+		0);
+	} else {
+	    bu_log("null nmg_vlblock_anim_upcall, no animation\n");
+	}
+	bn_vlblock_free(vbp);
+    }
+}
+
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

src/libnmg/nmg_index.c

@@ -0,0 +1,949 @@
+/*                     N M G _ I N D E X . C
+ * BRL-CAD
+ *
+ * Copyright (c) 1990-2016 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @addtogroup nmg */
+/** @{ */
+/** @file primitives/nmg/nmg_index.c
+ *
+ * Handle counting and re-indexing of NMG data structures.
+ *
+ */
+/** @} */
+
+#include "common.h"
+
+#include "nmg.h"
+
+#include "bio.h"
+#include "bu/log.h"
+#include "bu/malloc.h"
+
+#include <stddef.h>
+#include <string.h>
+#include <math.h>
+
+
+/**
+ * Return the structure index number of an arbitrary NMG structure.
+ *
+ * Returns -
+ * >=0 index number
+ * -1 pointed at struct bu_list embedded within NMG structure.
+ * -2 error:  unknown magic number
+ */
+int
+nmg_index_of_struct(register const uint32_t *p)
+{
+    switch (*p) {
+	case NMG_MODEL_MAGIC:
+	    return ((struct model *)p)->index;
+	case NMG_REGION_MAGIC:
+	    return ((struct nmgregion *)p)->index;
+	case NMG_REGION_A_MAGIC:
+	    return ((struct nmgregion_a *)p)->index;
+	case NMG_SHELL_MAGIC:
+	    return ((struct shell *)p)->index;
+	case NMG_SHELL_A_MAGIC:
+	    return ((struct shell_a *)p)->index;
+	case NMG_FACEUSE_MAGIC:
+	    return ((struct faceuse *)p)->index;
+	case NMG_FACE_MAGIC:
+	    return ((struct face *)p)->index;
+	case NMG_FACE_G_PLANE_MAGIC:
+	    return ((struct face_g_plane *)p)->index;
+	case NMG_FACE_G_SNURB_MAGIC:
+	    return ((struct face_g_snurb *)p)->index;
+	case NMG_LOOPUSE_MAGIC:
+	    return ((struct loopuse *)p)->index;
+	case NMG_LOOP_MAGIC:
+	    return ((struct loop *)p)->index;
+	case NMG_LOOP_G_MAGIC:
+	    return ((struct loop_g *)p)->index;
+	case NMG_EDGEUSE_MAGIC:
+	    return ((struct edgeuse *)p)->index;
+	case NMG_EDGEUSE2_MAGIC:
+	    /* Points to l2 inside edgeuse */
+	    return BU_LIST_MAIN_PTR(edgeuse, p, l2)->index;
+	case NMG_EDGE_MAGIC:
+	    return ((struct edge *)p)->index;
+	case NMG_EDGE_G_LSEG_MAGIC:
+	    return ((struct edge_g_lseg *)p)->index;
+	case NMG_EDGE_G_CNURB_MAGIC:
+	    return ((struct edge_g_cnurb *)p)->index;
+	case NMG_VERTEXUSE_MAGIC:
+	    return ((struct vertexuse *)p)->index;
+	case NMG_VERTEXUSE_A_PLANE_MAGIC:
+	    return ((struct vertexuse_a_plane *)p)->index;
+	case NMG_VERTEXUSE_A_CNURB_MAGIC:
+	    return ((struct vertexuse_a_cnurb *)p)->index;
+	case NMG_VERTEX_MAGIC:
+	    return ((struct vertex *)p)->index;
+	case NMG_VERTEX_G_MAGIC:
+	    return ((struct vertex_g *)p)->index;
+	case BU_LIST_HEAD_MAGIC:
+	    /* indicate special list head encountered */
+	    return -1;
+    }
+    /* default */
+    bu_log("nmg_index_of_struct: magicp = %p, magic = %x\n", (void *)p, *p);
+    return -2;	/* indicate error */
+}
+
+
+#define NMG_HIGH_BIT 0x80000000
+
+#define NMG_MARK_INDEX(_p)	((_p)->index |= NMG_HIGH_BIT)
+
+#define NMG_ASSIGN_NEW_INDEX(_p)	\
+    { if (((_p)->index & NMG_HIGH_BIT) != 0) \
+	    (_p)->index = newindex++; }
+
+/**
+ * Helper routine
+ */
+static void
+nmg_mark_edge_g(uint32_t *magic_p)
+{
+    if (!magic_p) bu_bomb("nmg_mark_edge_g bad magic\n");
+    switch (*magic_p) {
+	case NMG_EDGE_G_LSEG_MAGIC: {
+	    struct edge_g_lseg *lseg = (struct edge_g_lseg *)magic_p;
+	    NMG_MARK_INDEX(lseg);
+	    return;
+	}
+	case NMG_EDGE_G_CNURB_MAGIC: {
+	    struct edge_g_cnurb *cnurb = (struct edge_g_cnurb *)magic_p;
+	    NMG_MARK_INDEX(cnurb);
+	    return;
+	}
+    }
+    bu_bomb("nmg_mark_edge_g() unknown magic\n");
+}
+
+
+/**
+ * First pass:  just set the high bit on all index words
+ *
+ * This is a separate function largely for the benefit of global optimizers,
+ * which tended to blow their brains out on such a large subroutine.
+ */
+void
+nmg_m_set_high_bit(struct model *m)
+{
+    struct nmgregion *r;
+    struct shell *s;
+    struct faceuse *fu;
+    struct face *f;
+    struct loopuse *lu;
+    struct loop *l;
+    register struct edgeuse *eu;
+    struct edge *e;
+    register struct vertexuse *vu;
+
+#define MARK_VU(_vu) { \
+	struct vertex *v; \
+	NMG_CK_VERTEXUSE(_vu); \
+	NMG_MARK_INDEX(_vu); \
+	if (_vu->a.magic_p) switch (*_vu->a.magic_p) { \
+		case NMG_VERTEXUSE_A_PLANE_MAGIC: \
+		    NMG_MARK_INDEX(_vu->a.plane_p); \
+		    break; \
+		case NMG_VERTEXUSE_A_CNURB_MAGIC: \
+		    NMG_MARK_INDEX(_vu->a.cnurb_p); \
+		    break; \
+	    } \
+	v = _vu->v_p; \
+	NMG_CK_VERTEX(v); \
+	NMG_MARK_INDEX(v); \
+	if (v->vg_p) { \
+	    NMG_CK_VERTEX_G(v->vg_p); \
+	    NMG_MARK_INDEX(v->vg_p); \
+	} \
+    }
+
+    NMG_CK_MODEL(m);
+    NMG_MARK_INDEX(m);
+
+    for (BU_LIST_FOR(r, nmgregion, &m->r_hd)) {
+	NMG_CK_REGION(r);
+	NMG_MARK_INDEX(r);
+	if (r->ra_p) {
+	    NMG_CK_REGION_A(r->ra_p);
+	    NMG_MARK_INDEX(r->ra_p);
+	}
+	for (BU_LIST_FOR(s, shell, &r->s_hd)) {
+	    NMG_CK_SHELL(s);
+	    NMG_MARK_INDEX(s);
+	    if (s->sa_p) {
+		NMG_CK_SHELL_A(s->sa_p);
+		NMG_MARK_INDEX(s->sa_p);
+	    }
+	    /* Faces in shell */
+	    for (BU_LIST_FOR(fu, faceuse, &s->fu_hd)) {
+		NMG_CK_FACEUSE(fu);
+		NMG_MARK_INDEX(fu);
+		f = fu->f_p;
+		NMG_CK_FACE(f);
+		NMG_MARK_INDEX(f);
+		if (f->g.magic_p) switch (*f->g.magic_p) {
+			case NMG_FACE_G_PLANE_MAGIC:
+			    NMG_MARK_INDEX(f->g.plane_p);
+			    break;
+			case NMG_FACE_G_SNURB_MAGIC:
+			    NMG_MARK_INDEX(f->g.snurb_p);
+			    break;
+		    }
+		/* Loops in face */
+		for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
+		    NMG_CK_LOOPUSE(lu);
+		    NMG_MARK_INDEX(lu);
+		    l = lu->l_p;
+		    NMG_CK_LOOP(l);
+		    NMG_MARK_INDEX(l);
+		    if (l->lg_p) {
+			NMG_CK_LOOP_G(l->lg_p);
+			NMG_MARK_INDEX(l->lg_p);
+		    }
+		    if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+			/* Loop of Lone vertex */
+			vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
+			MARK_VU(vu);
+			continue;
+		    }
+		    for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+			NMG_CK_EDGEUSE(eu);
+			NMG_MARK_INDEX(eu);
+			e = eu->e_p;
+			NMG_CK_EDGE(e);
+			NMG_MARK_INDEX(e);
+			if (eu->g.magic_p) nmg_mark_edge_g(eu->g.magic_p);
+			vu = eu->vu_p;
+			MARK_VU(vu);
+		    }
+		}
+	    }
+	    /* Wire loops in shell */
+	    for (BU_LIST_FOR(lu, loopuse, &s->lu_hd)) {
+		NMG_CK_LOOPUSE(lu);
+		NMG_MARK_INDEX(lu);
+		l = lu->l_p;
+		NMG_CK_LOOP(l);
+		NMG_MARK_INDEX(l);
+		if (l->lg_p) {
+		    NMG_CK_LOOP_G(l->lg_p);
+		    NMG_MARK_INDEX(l->lg_p);
+		}
+		if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+		    /* Wire loop of Lone vertex */
+		    vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
+		    MARK_VU(vu);
+		    continue;
+		}
+		for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+		    NMG_CK_EDGEUSE(eu);
+		    NMG_MARK_INDEX(eu);
+		    e = eu->e_p;
+		    NMG_CK_EDGE(e);
+		    NMG_MARK_INDEX(e);
+		    if (eu->g.magic_p) nmg_mark_edge_g(eu->g.magic_p);
+		    vu = eu->vu_p;
+		    MARK_VU(vu);
+		}
+	    }
+	    /* Wire edges in shell */
+	    for (BU_LIST_FOR(eu, edgeuse, &s->eu_hd)) {
+		NMG_CK_EDGEUSE(eu);
+		NMG_MARK_INDEX(eu);
+		e = eu->e_p;
+		NMG_CK_EDGE(e);
+		NMG_MARK_INDEX(e);
+		if (eu->g.magic_p) nmg_mark_edge_g(eu->g.magic_p);
+		vu = eu->vu_p;
+		MARK_VU(vu);
+	    }
+	    /* Lone vertex in shell */
+	    vu = s->vu_p;
+	    if (vu) {
+		MARK_VU(vu);
+	    }
+	}
+    }
+#undef MARK_VU
+}
+
+
+/**
+ * Reassign index numbers to all the data structures in a model.
+ * The model structure will get index 0, all others will be sequentially
+ * assigned after that.
+ *
+ * Because the first pass has done extensive error checking,
+ * the second pass can do less.
+ */
+void
+nmg_m_reindex(struct model *m, register long int newindex)
+{
+    struct nmgregion *r;
+    struct shell *s;
+    struct faceuse *fu;
+    struct face *f;
+    struct loopuse *lu;
+    struct loop *l;
+    register struct edgeuse *eu;
+    struct edge *e;
+    register struct vertexuse *vu;
+    struct vertex *v;
+
+#define ASSIGN_VU(_vu) { \
+	NMG_CK_VERTEXUSE(_vu); \
+	NMG_ASSIGN_NEW_INDEX(_vu); \
+	if (_vu->a.magic_p) switch (*_vu->a.magic_p) { \
+		case NMG_VERTEXUSE_A_PLANE_MAGIC: \
+		    NMG_ASSIGN_NEW_INDEX(_vu->a.plane_p); \
+		    break; \
+		case NMG_VERTEXUSE_A_CNURB_MAGIC: \
+		    NMG_ASSIGN_NEW_INDEX(_vu->a.cnurb_p); \
+		    break; \
+	    } \
+	v = _vu->v_p; \
+	NMG_CK_VERTEX(v); \
+	NMG_ASSIGN_NEW_INDEX(v); \
+	if (v->vg_p) NMG_ASSIGN_NEW_INDEX(v->vg_p); \
+    }
+
+    NMG_CK_MODEL(m);
+    if (m->index != 0) bu_log("nmg_m_reindex() m->index=%ld\n", m->index);
+    if (newindex < 0) bu_log("nmg_m_reindex() newindex(%ld) < 0\n", newindex);
+
+    /* First pass:  set high bits */
+    nmg_m_set_high_bit(m);
+
+    /*
+     * Second pass:  assign new index number
+     */
+
+    NMG_ASSIGN_NEW_INDEX(m);
+    for (BU_LIST_FOR(r, nmgregion, &m->r_hd)) {
+	NMG_CK_REGION(r);
+	NMG_ASSIGN_NEW_INDEX(r);
+	if (r->ra_p) NMG_ASSIGN_NEW_INDEX(r->ra_p);
+	for (BU_LIST_FOR(s, shell, &r->s_hd)) {
+	    NMG_CK_SHELL(s);
+	    NMG_ASSIGN_NEW_INDEX(s);
+	    if (s->sa_p) NMG_ASSIGN_NEW_INDEX(s->sa_p);
+	    /* Faces in shell */
+	    for (BU_LIST_FOR(fu, faceuse, &s->fu_hd)) {
+		NMG_CK_FACEUSE(fu);
+		NMG_ASSIGN_NEW_INDEX(fu);
+		f = fu->f_p;
+		NMG_CK_FACE(f);
+		NMG_ASSIGN_NEW_INDEX(f);
+		if (f->g.plane_p) switch (*f->g.magic_p) {
+			case NMG_FACE_G_PLANE_MAGIC:
+			    NMG_ASSIGN_NEW_INDEX(f->g.plane_p);
+			    break;
+			case NMG_FACE_G_SNURB_MAGIC:
+			    NMG_ASSIGN_NEW_INDEX(f->g.snurb_p);
+			    break;
+		    }
+		/* Loops in face */
+		for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
+		    NMG_CK_LOOPUSE(lu);
+		    NMG_ASSIGN_NEW_INDEX(lu);
+		    l = lu->l_p;
+		    NMG_CK_LOOP(l);
+		    NMG_ASSIGN_NEW_INDEX(l);
+		    if (l->lg_p) NMG_ASSIGN_NEW_INDEX(l->lg_p);
+		    if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+			/* Loop of Lone vertex */
+			vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
+			ASSIGN_VU(vu);
+			continue;
+		    }
+		    for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+			NMG_CK_EDGEUSE(eu);
+			NMG_ASSIGN_NEW_INDEX(eu);
+			e = eu->e_p;
+			NMG_CK_EDGE(e);
+			NMG_ASSIGN_NEW_INDEX(e);
+			if (eu->g.magic_p) switch (*eu->g.magic_p) {
+				case NMG_EDGE_G_LSEG_MAGIC:
+				    NMG_ASSIGN_NEW_INDEX(eu->g.lseg_p);
+				    break;
+				case NMG_EDGE_G_CNURB_MAGIC:
+				    NMG_ASSIGN_NEW_INDEX(eu->g.cnurb_p);
+				    break;
+			    }
+			vu = eu->vu_p;
+			ASSIGN_VU(vu);
+		    }
+		}
+	    }
+	    /* Wire loops in shell */
+	    for (BU_LIST_FOR(lu, loopuse, &s->lu_hd)) {
+		NMG_CK_LOOPUSE(lu);
+		NMG_ASSIGN_NEW_INDEX(lu);
+		l = lu->l_p;
+		NMG_CK_LOOP(l);
+		NMG_ASSIGN_NEW_INDEX(l);
+		if (l->lg_p) NMG_ASSIGN_NEW_INDEX(l->lg_p);
+		if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+		    /* Wire loop of Lone vertex */
+		    vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
+		    ASSIGN_VU(vu);
+		    continue;
+		}
+		for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+		    NMG_CK_EDGEUSE(eu);
+		    NMG_ASSIGN_NEW_INDEX(eu);
+		    e = eu->e_p;
+		    NMG_CK_EDGE(e);
+		    NMG_ASSIGN_NEW_INDEX(e);
+		    if (eu->g.magic_p) switch (*eu->g.magic_p) {
+			    case NMG_EDGE_G_LSEG_MAGIC:
+				NMG_ASSIGN_NEW_INDEX(eu->g.lseg_p);
+				break;
+			    case NMG_EDGE_G_CNURB_MAGIC:
+				NMG_ASSIGN_NEW_INDEX(eu->g.cnurb_p);
+				break;
+			}
+		    vu = eu->vu_p;
+		    ASSIGN_VU(vu);
+		}
+	    }
+	    /* Wire edges in shell */
+	    for (BU_LIST_FOR(eu, edgeuse, &s->eu_hd)) {
+		NMG_CK_EDGEUSE(eu);
+		NMG_ASSIGN_NEW_INDEX(eu);
+		e = eu->e_p;
+		NMG_CK_EDGE(e);
+		NMG_ASSIGN_NEW_INDEX(e);
+		if (eu->g.magic_p) switch (*eu->g.magic_p) {
+			case NMG_EDGE_G_LSEG_MAGIC:
+			    NMG_ASSIGN_NEW_INDEX(eu->g.lseg_p);
+			    break;
+			case NMG_EDGE_G_CNURB_MAGIC:
+			    NMG_ASSIGN_NEW_INDEX(eu->g.cnurb_p);
+			    break;
+		    }
+		vu = eu->vu_p;
+		ASSIGN_VU(vu);
+	    }
+	    /* Lone vertex in shell */
+	    vu = s->vu_p;
+	    if (vu) {
+		ASSIGN_VU(vu);
+	    }
+	}
+    }
+#undef ASSIGN_VU
+
+    if (nmg_debug & DEBUG_BASIC) {
+	bu_log("nmg_m_reindex() oldmax=%ld, new%ld=>%ld\n",
+	       m->maxindex, m->index, newindex);
+    }
+    m->maxindex = newindex;
+}
+
+
+void
+nmg_vls_struct_counts(struct bu_vls *str, const struct nmg_struct_counts *ctr)
+{
+    BU_CK_VLS(str);
+
+    bu_vls_printf(str, " Actual structure counts:\n");
+    bu_vls_printf(str, "\t%6ld model\n", ctr->model);
+    bu_vls_printf(str, "\t%6ld region\n", ctr->region);
+    bu_vls_printf(str, "\t%6ld region_a\n", ctr->region_a);
+    bu_vls_printf(str, "\t%6ld shell\n", ctr->shell);
+    bu_vls_printf(str, "\t%6ld shell_a\n", ctr->shell_a);
+    bu_vls_printf(str, "\t%6ld face\n", ctr->face);
+    bu_vls_printf(str, "\t%6ld face_g_plane\n", ctr->face_g_plane);
+    bu_vls_printf(str, "\t%6ld face_g_snurb\n", ctr->face_g_snurb);
+    bu_vls_printf(str, "\t%6ld faceuse\n", ctr->faceuse);
+    bu_vls_printf(str, "\t%6ld loopuse\n", ctr->loopuse);
+    bu_vls_printf(str, "\t%6ld loop\n", ctr->loop);
+    bu_vls_printf(str, "\t%6ld loop_g\n", ctr->loop_g);
+    bu_vls_printf(str, "\t%6ld edgeuse\n", ctr->edgeuse);
+    bu_vls_printf(str, "\t%6ld edge\n", ctr->edge);
+    bu_vls_printf(str, "\t%6ld edge_g_lseg\n", ctr->edge_g_lseg);
+    bu_vls_printf(str, "\t%6ld edge_g_cnurb\n", ctr->edge_g_cnurb);
+    bu_vls_printf(str, "\t%6ld vertexuse\n", ctr->vertexuse);
+    bu_vls_printf(str, "\t%6ld vertexuse_a_plane\n", ctr->vertexuse_a_plane);
+    bu_vls_printf(str, "\t%6ld vertexuse_a_cnurb\n", ctr->vertexuse_a_cnurb);
+    bu_vls_printf(str, "\t%6ld vertex\n", ctr->vertex);
+    bu_vls_printf(str, "\t%6ld vertex_g\n", ctr->vertex_g);
+    bu_vls_printf(str, " Abstractions:\n");
+    bu_vls_printf(str, "\t%6ld max_structs\n", ctr->max_structs);
+    bu_vls_printf(str, "\t%6ld face_loops\n", ctr->face_loops);
+    bu_vls_printf(str, "\t%6ld face_edges\n", ctr->face_edges);
+    bu_vls_printf(str, "\t%6ld face_lone_verts\n", ctr->face_lone_verts);
+    bu_vls_printf(str, "\t%6ld wire_loops\n", ctr->wire_loops);
+    bu_vls_printf(str, "\t%6ld wire_loop_edges\n", ctr->wire_loop_edges);
+    bu_vls_printf(str, "\t%6ld wire_edges\n", ctr->wire_edges);
+    bu_vls_printf(str, "\t%6ld wire_lone_verts\n", ctr->wire_lone_verts);
+    bu_vls_printf(str, "\t%6ld shells_of_lone_vert\n", ctr->shells_of_lone_vert);
+}
+
+
+void
+nmg_pr_struct_counts(const struct nmg_struct_counts *ctr, const char *str)
+{
+    struct bu_vls vls = BU_VLS_INIT_ZERO;
+
+    bu_log("nmg_pr_count(%s)\n", str);
+
+    nmg_vls_struct_counts(&vls, ctr);
+    bu_log("%s", bu_vls_addr(&vls));
+    bu_vls_free(&vls);
+}
+
+
+/**
+ * Returns -
+ * Pointer to magic-number/structure-base pointer array,
+ * indexed by nmg structure index.
+ * Caller is responsible for freeing it.
+ */
+uint32_t **
+nmg_m_struct_count(register struct nmg_struct_counts *ctr, const struct model *m)
+{
+    struct nmgregion *r;
+    struct shell *s;
+    struct faceuse *fu;
+    struct face *f;
+    struct loopuse *lu;
+    struct loop *l;
+    struct edgeuse *eu;
+    struct edge *e;
+    struct vertexuse *vu;
+    struct vertex *v;
+    register uint32_t **ptrs;
+
+#define NMG_UNIQ_INDEX(_p, _type)	\
+    if ((_p)->index > m->maxindex) { \
+	bu_log("%p (%s) has index %ld, m->maxindex=%ld\n", (void *)(_p), \
+	       bu_identify_magic(*((uint32_t *)(_p))), (_p)->index, m->maxindex); \
+	bu_bomb("nmg_m_struct_count index overflow\n"); \
+    } \
+    if (ptrs[(_p)->index] == NULL) { \
+	ptrs[(_p)->index] = (uint32_t *)(_p); \
+	ctr->_type++; \
+    }
+
+#define UNIQ_VU(_vu) { \
+	NMG_CK_VERTEXUSE(_vu); \
+	NMG_UNIQ_INDEX(_vu, vertexuse); \
+	if (_vu->a.magic_p) switch (*_vu->a.magic_p) { \
+		case NMG_VERTEXUSE_A_PLANE_MAGIC: \
+		    NMG_UNIQ_INDEX(_vu->a.plane_p, vertexuse_a_plane); \
+		    break; \
+		case NMG_VERTEXUSE_A_CNURB_MAGIC: \
+		    NMG_UNIQ_INDEX(_vu->a.cnurb_p, vertexuse_a_cnurb); \
+		    break; \
+	    } \
+	v = _vu->v_p; \
+	NMG_CK_VERTEX(v); \
+	NMG_UNIQ_INDEX(v, vertex); \
+	if (v->vg_p) { \
+	    NMG_CK_VERTEX_G(v->vg_p); \
+	    NMG_UNIQ_INDEX(v->vg_p, vertex_g); \
+	} \
+    }
+
+    NMG_CK_MODEL(m);
+    memset((char *)ctr, 0, sizeof(*ctr));
+
+    ptrs = (uint32_t **)bu_calloc(m->maxindex+1, sizeof(uint32_t *), "nmg_m_count ptrs[]");
+
+    NMG_UNIQ_INDEX(m, model);
+    ctr->max_structs = m->maxindex;
+    for (BU_LIST_FOR(r, nmgregion, &m->r_hd)) {
+	NMG_CK_REGION(r);
+	NMG_UNIQ_INDEX(r, region);
+	if (r->ra_p) {
+	    NMG_CK_REGION_A(r->ra_p);
+	    NMG_UNIQ_INDEX(r->ra_p, region_a);
+	}
+	for (BU_LIST_FOR(s, shell, &r->s_hd)) {
+	    NMG_CK_SHELL(s);
+	    NMG_UNIQ_INDEX(s, shell);
+	    if (s->sa_p) {
+		NMG_CK_SHELL_A(s->sa_p);
+		NMG_UNIQ_INDEX(s->sa_p, shell_a);
+	    }
+	    /* Faces in shell */
+	    for (BU_LIST_FOR(fu, faceuse, &s->fu_hd)) {
+		NMG_CK_FACEUSE(fu);
+		NMG_UNIQ_INDEX(fu, faceuse);
+		f = fu->f_p;
+		NMG_CK_FACE(f);
+		NMG_UNIQ_INDEX(f, face);
+		if (f->g.magic_p) switch (*f->g.magic_p) {
+			case NMG_FACE_G_PLANE_MAGIC:
+			    NMG_UNIQ_INDEX(f->g.plane_p, face_g_plane);
+			    break;
+			case NMG_FACE_G_SNURB_MAGIC:
+			    NMG_UNIQ_INDEX(f->g.snurb_p, face_g_snurb);
+			    break;
+		    }
+		/* Loops in face */
+		for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
+		    NMG_CK_LOOPUSE(lu);
+		    NMG_UNIQ_INDEX(lu, loopuse);
+		    l = lu->l_p;
+		    NMG_CK_LOOP(l);
+		    NMG_UNIQ_INDEX(l, loop);
+		    if (l->lg_p) {
+			NMG_CK_LOOP_G(l->lg_p);
+			NMG_UNIQ_INDEX(l->lg_p, loop_g);
+		    }
+		    if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+			/* Loop of Lone vertex */
+			ctr->face_lone_verts++;
+			vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
+			UNIQ_VU(vu);
+			continue;
+		    }
+		    ctr->face_loops++;
+		    for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+			ctr->face_edges++;
+			NMG_CK_EDGEUSE(eu);
+			NMG_UNIQ_INDEX(eu, edgeuse);
+			e = eu->e_p;
+			NMG_CK_EDGE(e);
+			NMG_UNIQ_INDEX(e, edge);
+			if (eu->g.magic_p) switch (*eu->g.magic_p) {
+				case NMG_EDGE_G_LSEG_MAGIC:
+				    NMG_UNIQ_INDEX(eu->g.lseg_p, edge_g_lseg);
+				    break;
+				case NMG_EDGE_G_CNURB_MAGIC:
+				    NMG_UNIQ_INDEX(eu->g.cnurb_p, edge_g_cnurb);
+				    break;
+			    }
+			vu = eu->vu_p;
+			UNIQ_VU(vu);
+		    }
+		}
+	    }
+	    /* Wire loops in shell */
+	    for (BU_LIST_FOR(lu, loopuse, &s->lu_hd)) {
+		NMG_CK_LOOPUSE(lu);
+		NMG_UNIQ_INDEX(lu, loopuse);
+		l = lu->l_p;
+		NMG_CK_LOOP(l);
+		NMG_UNIQ_INDEX(l, loop);
+		if (l->lg_p) {
+		    NMG_CK_LOOP_G(l->lg_p);
+		    NMG_UNIQ_INDEX(l->lg_p, loop_g);
+		}
+		if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+		    ctr->wire_lone_verts++;
+		    /* Wire loop of Lone vertex */
+		    vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
+		    UNIQ_VU(vu);
+		    continue;
+		}
+		ctr->wire_loops++;
+		for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+		    NMG_CK_EDGEUSE(eu);
+		    NMG_UNIQ_INDEX(eu, edgeuse);
+		    e = eu->e_p;
+		    NMG_CK_EDGE(e);
+		    NMG_UNIQ_INDEX(e, edge);
+		    if (eu->g.magic_p) switch (*eu->g.magic_p) {
+			    case NMG_EDGE_G_LSEG_MAGIC:
+				NMG_UNIQ_INDEX(eu->g.lseg_p, edge_g_lseg);
+				break;
+			    case NMG_EDGE_G_CNURB_MAGIC:
+				NMG_UNIQ_INDEX(eu->g.cnurb_p, edge_g_cnurb);
+				break;
+			}
+		    vu = eu->vu_p;
+		    UNIQ_VU(vu);
+		    ctr->wire_loop_edges++;
+		}
+	    }
+	    /* Wire edges in shell */
+	    for (BU_LIST_FOR(eu, edgeuse, &s->eu_hd)) {
+		NMG_CK_EDGEUSE(eu);
+		ctr->wire_edges++;
+		NMG_UNIQ_INDEX(eu, edgeuse);
+		e = eu->e_p;
+		NMG_CK_EDGE(e);
+		NMG_UNIQ_INDEX(e, edge);
+		if (eu->g.magic_p) switch (*eu->g.magic_p) {
+			case NMG_EDGE_G_LSEG_MAGIC:
+			    NMG_UNIQ_INDEX(eu->g.lseg_p, edge_g_lseg);
+			    break;
+			case NMG_EDGE_G_CNURB_MAGIC:
+			    NMG_UNIQ_INDEX(eu->g.cnurb_p, edge_g_cnurb);
+			    break;
+		    }
+		vu = eu->vu_p;
+		UNIQ_VU(vu);
+	    }
+	    /* Lone vertex in shell */
+	    vu = s->vu_p;
+	    if (vu) {
+		ctr->shells_of_lone_vert++;
+		UNIQ_VU(vu);
+	    }
+	}
+    }
+    /* Caller must free them */
+    return ptrs;
+#undef UNIQ_VU
+}
+
+
+/**
+ *
+ * Count 'em up, and print 'em out.
+ */
+void
+nmg_pr_m_struct_counts(const struct model *m, const char *str)
+{
+    struct nmg_struct_counts cnts;
+    uint32_t **tab;
+
+    NMG_CK_MODEL(m);
+
+    tab = nmg_m_struct_count(&cnts, m);
+    bu_free((void *)tab, "nmg_m_struct_count");
+    nmg_pr_struct_counts(&cnts, str);
+}
+
+
+/**
+ * Combine two NMG model trees into one single NMG model.  The
+ * first model inherits the nmgregions of the second.  The second
+ * model pointer is freed before return.
+ */
+void
+nmg_merge_models(struct model *m1, struct model *m2)
+{
+    struct nmgregion *r;
+
+    NMG_CK_MODEL(m1);
+    NMG_CK_MODEL(m2);
+
+    /* first reorder the first model to "compress" the
+     * number space if possible.
+     */
+    nmg_m_reindex(m1, 0);
+
+    if (m1 == m2) /* nothing to do */
+	return;
+
+    /* now re-order the second model starting with an index number
+     * of m1->maxindex.
+     *
+     * We might get away with using m1->maxindex-1, since the first
+     * value is assigned to the second model structure, and we will
+     * shortly be freeing the second model struct.
+     */
+
+    nmg_m_reindex(m2, m1->maxindex);
+    m1->maxindex = m2->maxindex;		/* big enough for both */
+
+    /* Rehome all the regions in m2, and move them from m2 to m1 */
+    for (BU_LIST_FOR(r, nmgregion, &(m2->r_hd))) {
+	NMG_CK_REGION(r);
+	r->m_p = m1;
+    }
+    BU_LIST_APPEND_LIST(&(m1->r_hd), &(m2->r_hd));
+
+    /* If there are any manifold tables, when the models are
+     * merged they become invalid and need to be regenerated.
+     * To avoid confusion, free them here so it is known they
+     * need to be regenerated.
+     */
+    if (m1->manifolds) {
+	bu_free((char *)m1->manifolds, "free manifolds table");
+	m1->manifolds = (char *)NULL;
+    }
+    if (m2->manifolds) {
+	bu_free((char *)m2->manifolds, "free manifolds table");
+	m2->manifolds = (char *)NULL;
+    }
+
+    FREE_MODEL(m2);
+}
+
+
+#define CHECK_INDEX(_p) if ((_p)->index > maxindex) maxindex = (_p)->index
+#define CHECK_VU_INDEX(_vu) {\
+	NMG_CK_VERTEXUSE(_vu); \
+	CHECK_INDEX(_vu); \
+	if (_vu->a.magic_p) switch (*_vu->a.magic_p) { \
+		case NMG_VERTEXUSE_A_PLANE_MAGIC: \
+		    CHECK_INDEX(_vu->a.plane_p); \
+		    break; \
+		case NMG_VERTEXUSE_A_CNURB_MAGIC: \
+		    CHECK_INDEX(_vu->a.cnurb_p); \
+		    break; \
+	    } \
+	v = _vu->v_p; \
+	NMG_CK_VERTEX(v); \
+	CHECK_INDEX(v); \
+	if (v->vg_p) CHECK_INDEX(v->vg_p); \
+    }
+
+long
+nmg_find_max_index(const struct model *m)
+{
+    long maxindex=0;
+    struct nmgregion *r;
+    struct shell *s;
+    struct faceuse *fu;
+    struct face *f;
+    struct loopuse *lu;
+    struct loop *l;
+    register struct edgeuse *eu;
+    struct edge *e;
+    register struct vertexuse *vu;
+    struct vertex *v;
+
+    for (BU_LIST_FOR(r, nmgregion, &m->r_hd)) {
+	NMG_CK_REGION(r);
+	CHECK_INDEX(r);
+	if (r->ra_p) CHECK_INDEX(r->ra_p);
+	for (BU_LIST_FOR(s, shell, &r->s_hd)) {
+	    NMG_CK_SHELL(s);
+	    CHECK_INDEX(s);
+	    if (s->sa_p) CHECK_INDEX(s->sa_p);
+	    /* Faces in shell */
+	    for (BU_LIST_FOR(fu, faceuse, &s->fu_hd)) {
+		NMG_CK_FACEUSE(fu);
+		CHECK_INDEX(fu);
+		f = fu->f_p;
+		NMG_CK_FACE(f);
+		CHECK_INDEX(f);
+		if (f->g.plane_p) switch (*f->g.magic_p) {
+			case NMG_FACE_G_PLANE_MAGIC:
+			    CHECK_INDEX(f->g.plane_p);
+			    break;
+			case NMG_FACE_G_SNURB_MAGIC:
+			    CHECK_INDEX(f->g.snurb_p);
+			    break;
+		    }
+		/* Loops in face */
+		for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
+		    NMG_CK_LOOPUSE(lu);
+		    CHECK_INDEX(lu);
+		    l = lu->l_p;
+		    NMG_CK_LOOP(l);
+		    CHECK_INDEX(l);
+		    if (l->lg_p) CHECK_INDEX(l->lg_p);
+		    if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+			/* Loop of Lone vertex */
+			vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
+			CHECK_VU_INDEX(vu);
+			continue;
+		    }
+		    for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+			NMG_CK_EDGEUSE(eu);
+			CHECK_INDEX(eu);
+			e = eu->e_p;
+			NMG_CK_EDGE(e);
+			CHECK_INDEX(e);
+			if (eu->g.magic_p) switch (*eu->g.magic_p) {
+				case NMG_EDGE_G_LSEG_MAGIC:
+				    CHECK_INDEX(eu->g.lseg_p);
+				    break;
+				case NMG_EDGE_G_CNURB_MAGIC:
+				    CHECK_INDEX(eu->g.cnurb_p);
+				    break;
+			    }
+			vu = eu->vu_p;
+			CHECK_VU_INDEX(vu);
+		    }
+		}
+	    }
+	    /* Wire loops in shell */
+	    for (BU_LIST_FOR(lu, loopuse, &s->lu_hd)) {
+		NMG_CK_LOOPUSE(lu);
+		CHECK_INDEX(lu);
+		l = lu->l_p;
+		NMG_CK_LOOP(l);
+		CHECK_INDEX(l);
+		if (l->lg_p) CHECK_INDEX(l->lg_p);
+		if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+		    /* Wire loop of Lone vertex */
+		    vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
+		    CHECK_VU_INDEX(vu);
+		    continue;
+		}
+		for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+		    NMG_CK_EDGEUSE(eu);
+		    CHECK_INDEX(eu);
+		    e = eu->e_p;
+		    NMG_CK_EDGE(e);
+		    CHECK_INDEX(e);
+		    if (eu->g.magic_p) switch (*eu->g.magic_p) {
+			    case NMG_EDGE_G_LSEG_MAGIC:
+				CHECK_INDEX(eu->g.lseg_p);
+				break;
+			    case NMG_EDGE_G_CNURB_MAGIC:
+				CHECK_INDEX(eu->g.cnurb_p);
+				break;
+			}
+		    vu = eu->vu_p;
+		    CHECK_VU_INDEX(vu);
+		}
+	    }
+	    /* Wire edges in shell */
+	    for (BU_LIST_FOR(eu, edgeuse, &s->eu_hd)) {
+		NMG_CK_EDGEUSE(eu);
+		CHECK_INDEX(eu);
+		e = eu->e_p;
+		NMG_CK_EDGE(e);
+		CHECK_INDEX(e);
+		if (eu->g.magic_p) switch (*eu->g.magic_p) {
+			case NMG_EDGE_G_LSEG_MAGIC:
+			    CHECK_INDEX(eu->g.lseg_p);
+			    break;
+			case NMG_EDGE_G_CNURB_MAGIC:
+			    CHECK_INDEX(eu->g.cnurb_p);
+			    break;
+		    }
+		vu = eu->vu_p;
+		CHECK_VU_INDEX(vu);
+	    }
+	    /* Lone vertex in shell */
+	    vu = s->vu_p;
+	    if (vu) {
+		CHECK_VU_INDEX(vu);
+	    }
+	}
+    }
+    return maxindex;
+}
+
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

src/libnmg/nmg_manif.c

@@ -0,0 +1,461 @@
+/*                     N M G _ M A N I F . C
+ * BRL-CAD
+ *
+ * Copyright (c) 1994-2016 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @addtogroup nmg */
+/** @{ */
+/** @file primitives/nmg/nmg_manif.c
+ *
+ * Routines for assessing the manifold dimension of an object.
+ *
+ */
+/** @} */
+
+#include "common.h"
+
+#include "nmg.h"
+#include "nmg_func.h"
+
+#include "bio.h"
+#include "bn/tol.h"
+#include "bu/log.h"
+#include "bu/malloc.h"
+
+#include <string.h>
+
+
+#define PAINT_INTERIOR 1
+#define PAINT_EXTERIOR 0
+
+#define BU_LIST_LINK_CHECK(p) \
+    if (BU_LIST_PNEXT_PLAST(bu_list, p) != p || \
+	BU_LIST_PLAST_PNEXT(bu_list, p) != p) { \
+	bu_log("%s[%d]: linked list integrity check failed\n", \
+	       __FILE__, __LINE__); \
+	bu_log("%p->forw(%p)->back = %p\n", \
+	       (void *)(p), (void *)(p)->forw, (void *)(p)->forw->back); \
+	bu_log("%p->back(%p)->forw = %p\n", \
+	       (void *)(p), (void *)(p)->back, (void *)(p)->back->forw); \
+	bu_bomb("Goodbye\n"); \
+    }
+
+
+/**
+ * Determine if a face has any "dangling" edges.
+ *
+ * Return
+ * 1 face has dangling edge
+ * 0 face does not have a dangling edge
+ */
+int
+nmg_dangling_face(const struct faceuse *fu, register const char *manifolds)
+{
+    struct loopuse *lu;
+    struct edgeuse *eu;
+    const struct edgeuse *eur;
+    struct faceuse *newfu;
+
+    NMG_CK_FACEUSE(fu);
+
+    if (nmg_debug & DEBUG_MANIF)
+	bu_log("nmg_dangling_face(%p %s)\n", (void *)fu, manifolds);
+
+    for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
+	NMG_CK_LOOPUSE(lu);
+	BU_LIST_LINK_CHECK(&lu->l);
+
+	if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_EDGEUSE_MAGIC) {
+	    /* go looking around each edge for a face of the same
+	     * shell which isn't us and isn't our mate.  If we
+	     * find us or our mate before another face of this
+	     * shell, we are non-3-manifold.
+	     */
+
+	    for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+
+		NMG_CK_EDGEUSE(eu);
+		BU_LIST_LINK_CHECK(&eu->l);
+
+		eur = nmg_radial_face_edge_in_shell(eu);
+		newfu = eur->up.lu_p->up.fu_p;
+
+		/* skip any known dangling-edge faces or
+		 * faces known to be 2manifolds.
+		 */
+		while (manifolds &&
+		       NMG_MANIFOLDS(manifolds, newfu) & NMG_2MANIFOLD &&
+		       eur != eu->eumate_p) {
+		    eur = nmg_radial_face_edge_in_shell(
+			eur->eumate_p);
+		    newfu = eur->up.lu_p->up.fu_p;
+		}
+
+		if (eur == eu->eumate_p) {
+		    goto out;
+		}
+	    }
+	}
+    }
+    eur = (const struct edgeuse *)NULL;
+
+out:
+    if (nmg_debug & DEBUG_BASIC) {
+	struct bn_tol tol;	/* HACK */
+	tol.magic = BN_TOL_MAGIC;
+	tol.dist = 1;
+	tol.dist_sq = tol.dist * tol.dist;
+	tol.perp = 1e-5;
+	tol.para = 1 - tol.perp;
+
+	bu_log("nmg_dangling_face(fu=%p, manifolds=%s) dangling_eu=%p\n", (void *)fu, manifolds, (void *)eur);
+	if (eur) nmg_pr_fu_around_eu(eur, &tol);
+    }
+    if ((nmg_debug & DEBUG_MANIF) && (eur != (const struct edgeuse *)NULL))
+	bu_log("\tdangling eu %p\n", (void *)eur);
+
+    return eur != (const struct edgeuse *)NULL;
+}
+
+
+/**
+ * "Paint" the elements of a face with a meaning.  For example
+ * mark everything in a face as being part of a 2-manifold
+ */
+static void paint_face(struct faceuse *fu, long *paint_table, long paint_color, char *paint_meaning, char *tbl)
+{
+    struct faceuse *newfu;
+    struct loopuse *lu;
+    struct edgeuse *eu;
+    const struct edgeuse *eur;
+
+    if (nmg_debug & DEBUG_MANIF)
+	bu_log("nmg_paint_face(%p, %ld)\n", (void *)fu, paint_color);
+
+    if (NMG_INDEX_VALUE(paint_table, fu->index) != 0)
+	return;
+
+    NMG_INDEX_ASSIGN(paint_table, fu, paint_color);
+
+    for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
+
+	if (nmg_debug & DEBUG_MANIF)
+	    bu_log("\tlu=%p\n", (void *)lu);
+
+	NMG_CK_LOOPUSE(lu);
+	BU_LIST_LINK_CHECK(&lu->l);
+
+	if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
+	    continue;
+
+	for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+	    if (nmg_debug & DEBUG_MANIF)
+		bu_log("\t\teu=%p\n", (void *)eu);
+	    NMG_CK_EDGEUSE(eu);
+	    NMG_CK_EDGEUSE(eu->eumate_p);
+	    eur = nmg_radial_face_edge_in_shell(eu);
+	    NMG_CK_EDGEUSE(eur);
+	    NMG_CK_FACEUSE(eur->up.lu_p->up.fu_p);
+	    newfu = eur->up.lu_p->up.fu_p;
+
+	    if (nmg_debug & DEBUG_MANIF)
+		bu_log("\t\t\teur=%p, newfu=%p\n", (void *)eur, (void *)newfu);
+
+	    BU_LIST_LINK_CHECK(&eu->l);
+	    BU_LIST_LINK_CHECK(&eur->l);
+
+	    while (NMG_MANIFOLDS(tbl, newfu) & NMG_2MANIFOLD &&
+		   eur != eu->eumate_p) {
+		eur = nmg_radial_face_edge_in_shell(
+		    eur->eumate_p);
+		newfu = eur->up.lu_p->up.fu_p;
+
+		if (nmg_debug & DEBUG_MANIF)
+		    bu_log("\t\t\teur=%p, newfu=%p\n", (void *)eur, (void *)newfu);
+	    }
+
+	    if (newfu == fu->fumate_p)
+		continue;
+	    else if (newfu->orientation == OT_SAME)
+		paint_face(newfu, paint_table, paint_color,
+			   paint_meaning, tbl);
+	    else {
+		/* mark this group as being interior */
+		paint_meaning[paint_color] = PAINT_INTERIOR;
+
+		if (nmg_debug & DEBUG_MANIF)
+		    bu_log("\t---- Painting fu %p as interior, new_fu = %p, eu=%p, eur=%p\n",
+			   (void *)fu, (void *)newfu, (void *)eu, (void *)eur);
+	    }
+	}
+    }
+}
+
+
+static void set_edge_sub_manifold(char *tbl, struct edgeuse *eu_p, char manifold)
+{
+
+    NMG_CK_EDGEUSE(eu_p);
+    NMG_CK_EDGE(eu_p->e_p);
+    NMG_CK_VERTEXUSE(eu_p->vu_p);
+    NMG_CK_VERTEX(eu_p->vu_p->v_p);
+
+    NMG_SET_MANIFOLD(tbl, eu_p, manifold);
+    NMG_SET_MANIFOLD(tbl, eu_p->e_p, manifold);
+
+    NMG_SET_MANIFOLD(tbl, eu_p->vu_p, manifold);
+    NMG_SET_MANIFOLD(tbl, eu_p->vu_p->v_p, manifold);
+}
+
+
+static void set_loop_sub_manifold(char *tbl, struct loopuse *lu_p, char manifold)
+{
+    struct edgeuse *eu_p;
+    struct vertexuse *vu_p;
+
+    NMG_CK_LOOPUSE(lu_p);
+
+    NMG_SET_MANIFOLD(tbl, lu_p, manifold);
+    NMG_SET_MANIFOLD(tbl, lu_p->l_p, manifold);
+    if (BU_LIST_FIRST_MAGIC(&lu_p->down_hd) == NMG_VERTEXUSE_MAGIC) {
+	vu_p = BU_LIST_FIRST(vertexuse, &lu_p->down_hd);
+	NMG_SET_MANIFOLD(tbl, vu_p, manifold);
+	NMG_SET_MANIFOLD(tbl, vu_p->v_p, manifold);
+    } else if (BU_LIST_FIRST_MAGIC(&lu_p->down_hd) == NMG_EDGEUSE_MAGIC) {
+	for (BU_LIST_FOR(eu_p, edgeuse, &lu_p->down_hd)) {
+	    BU_LIST_LINK_CHECK(&eu_p->l);
+	    set_edge_sub_manifold(tbl, eu_p, manifold);
+	}
+    } else
+	bu_bomb("bad child of loopuse\n");
+}
+static void set_face_sub_manifold(char *tbl, struct faceuse *fu_p, char manifold)
+{
+    struct loopuse *lu_p;
+
+    NMG_CK_FACEUSE(fu_p);
+    NMG_CK_FACE(fu_p->f_p);
+
+    NMG_SET_MANIFOLD(tbl, fu_p, manifold);
+    NMG_SET_MANIFOLD(tbl, fu_p->f_p, manifold);
+    for (BU_LIST_FOR(lu_p, loopuse, &fu_p->lu_hd)) {
+	BU_LIST_LINK_CHECK(&lu_p->l);
+	set_loop_sub_manifold(tbl, lu_p, manifold);
+    }
+}
+
+
+char *
+nmg_shell_manifolds(struct shell *sp, char *tbl)
+{
+    struct edgeuse *eu_p;
+    struct loopuse *lu_p;
+    struct faceuse *fu_p;
+    char *paint_meaning;
+    long *paint_table;
+    long paint_color;
+    int found;
+
+    if (nmg_debug & DEBUG_MANIF)
+	bu_log("nmg_shell_manifolds(%p)\n", (void *)sp);
+
+    NMG_CK_SHELL(sp);
+
+    if (tbl == (char *)NULL)
+	tbl = (char *)bu_calloc(sp->r_p->m_p->maxindex, 1, "manifold table");
+
+    /*
+     * points in shells form 0-manifold objects.
+     */
+    if (sp->vu_p) {
+	NMG_CK_VERTEXUSE(sp->vu_p);
+	NMG_SET_MANIFOLD(tbl, sp, NMG_0MANIFOLD);
+	NMG_SET_MANIFOLD(tbl, sp->vu_p, NMG_0MANIFOLD);
+	NMG_SET_MANIFOLD(tbl, sp->vu_p->v_p, NMG_0MANIFOLD);
+	BU_LIST_LINK_CHECK(&sp->vu_p->l);
+    }
+
+    /* edges in shells are (components of)
+     * 1-manifold objects.
+     */
+    if (BU_LIST_NON_EMPTY(&sp->eu_hd)) {
+
+	for (BU_LIST_FOR(eu_p, edgeuse, &sp->eu_hd)) {
+	    BU_LIST_LINK_CHECK(&eu_p->l);
+	    set_edge_sub_manifold(tbl, eu_p, NMG_1MANIFOLD);
+	}
+	NMG_SET_MANIFOLD(tbl, sp, NMG_1MANIFOLD);
+    }
+
+    /* loops in shells are (components of)
+     * 1-manifold objects.
+     */
+    if (BU_LIST_NON_EMPTY(&sp->lu_hd)) {
+
+	for (BU_LIST_FOR(lu_p, loopuse, &sp->lu_hd)) {
+	    BU_LIST_LINK_CHECK(&lu_p->l);
+
+	    set_loop_sub_manifold(tbl, lu_p, NMG_1MANIFOLD);
+	}
+	NMG_SET_MANIFOLD(tbl, sp, NMG_1MANIFOLD);
+    }
+
+    if (nmg_debug & DEBUG_MANIF)
+	bu_log("starting manifold classification on shell faces\n");
+
+    /*
+     * faces may be either 2 or 3 manifold components.
+     */
+    if (BU_LIST_IS_EMPTY(&sp->fu_hd))
+	return tbl;
+
+
+    /* mark all externally dangling faces and faces
+     * with or adjacent to dangling edges.
+     */
+    do {
+	found = 0;
+	for (BU_LIST_FOR(fu_p, faceuse, &sp->fu_hd)) {
+	    NMG_CK_FACEUSE(fu_p);
+	    BU_LIST_LINK_CHECK(&fu_p->l);
+
+	    /* if this has already been marked as a 2-manifold
+	     * then we don't need to check it again
+	     */
+	    if (NMG_MANIFOLDS(tbl, fu_p) & NMG_2MANIFOLD)
+		continue;
+
+	    if (nmg_dangling_face(fu_p, tbl)) {
+		found = 1;
+
+		NMG_SET_MANIFOLD(tbl, fu_p, NMG_2MANIFOLD);
+
+		if (nmg_debug & DEBUG_MANIF)
+		    bu_log("found dangling face\n");
+	    }
+	}
+    } while (found);
+
+    /* paint the exterior faces to discover what the
+     * actual enclosing shell is
+     */
+
+    if (nmg_debug & DEBUG_MANIF)
+	bu_log("starting to paint non-dangling faces\n");
+
+    paint_meaning = (char *)bu_calloc(sp->r_p->m_p->maxindex, sizeof(char), "paint meaning table");
+    paint_table = (long *)bu_calloc(sp->r_p->m_p->maxindex, sizeof(long), "paint table");
+    paint_color = 1;
+
+    for (BU_LIST_FOR(fu_p, faceuse, &sp->fu_hd)) {
+	BU_LIST_LINK_CHECK(&fu_p->l);
+
+	if (fu_p->orientation != OT_SAME ||
+	    NMG_INDEX_VALUE(paint_table, fu_p->index) != 0)
+	    continue;
+
+	paint_face(fu_p, paint_table, paint_color, paint_meaning, tbl);
+	paint_color++;
+    }
+
+    if (nmg_debug & DEBUG_MANIF)
+	bu_log("painting done, looking at colors\n");
+
+    /* all the faces painted with "interior" paint are 2manifolds
+     * those faces still painted with "exterior" paint are
+     * 3manifolds, i.e. part of the enclosing surface
+     */
+    for (BU_LIST_FOR(fu_p, faceuse, &sp->fu_hd)) {
+	BU_LIST_LINK_CHECK(&fu_p->l);
+	paint_color = NMG_INDEX_VALUE((long *)paint_table, fu_p->index);
+
+	/* this should never trigger. */
+	/* it is easier to compare against the max model index */
+	if (paint_color > sp->r_p->m_p->maxindex - 1) {
+	    bu_log("nmg_shell_manifolds(): ERROR, color index out of range (%ld > %ld)\n",
+		   paint_color, sp->r_p->m_p->maxindex - 1);
+	    bu_bomb("nmg_shell_manifolds(): ERROR, color index out of range\n");
+	    break;
+	}
+
+	if (NMG_INDEX_VALUE(paint_meaning, paint_color) == PAINT_INTERIOR) {
+	    set_face_sub_manifold(tbl, fu_p, NMG_2MANIFOLD);
+	} else if (NMG_INDEX_VALUE(paint_meaning, paint_color) == PAINT_EXTERIOR) {
+	    set_face_sub_manifold(tbl, fu_p, NMG_3MANIFOLD);
+	}
+    }
+
+    bu_free(paint_meaning, "paint meaning table");
+    bu_free(paint_table, "paint table");
+
+    for (BU_LIST_FOR(fu_p, faceuse, &sp->fu_hd)) {
+	BU_LIST_LINK_CHECK(&fu_p->l);
+
+	if (fu_p->orientation != OT_SAME)
+	    NMG_CP_MANIFOLD(tbl, fu_p, fu_p->fumate_p);
+    }
+
+    return tbl;
+}
+
+
+char *
+nmg_manifolds(struct model *m)
+{
+    struct nmgregion *rp;
+    struct shell *sp;
+    char *tbl;
+
+
+    NMG_CK_MODEL(m);
+    if (nmg_debug & DEBUG_MANIF)
+	bu_log("nmg_manifolds(%p)\n", (void *)m);
+
+    tbl = (char *)bu_calloc(m->maxindex, 1, "manifold table");
+
+
+    for (BU_LIST_FOR(rp, nmgregion, &m->r_hd)) {
+	NMG_CK_REGION(rp);
+	BU_LIST_LINK_CHECK(&rp->l);
+
+	for (BU_LIST_FOR(sp, shell, &rp->s_hd)) {
+
+	    NMG_CK_SHELL(sp);
+	    BU_LIST_LINK_CHECK(&sp->l);
+
+	    nmg_shell_manifolds(sp, tbl);
+
+	    /* make sure the region manifold bits are a superset
+	     * of the shell manifold bits
+	     */
+	    NMG_CP_MANIFOLD(tbl, rp, sp);
+	}
+    }
+
+    return tbl;
+}
+
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

src/libnmg/nmg_mesh.c

@@ -0,0 +1,669 @@
+/*                      N M G _ M E S H . C
+ * BRL-CAD
+ *
+ * Copyright (c) 1989-2016 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @addtogroup nmg */
+/** @{ */
+/** @file primitives/nmg/nmg_mesh.c
+ *
+ * Meshing routines for n-Manifold Geometry.
+ * This stuff is destined to be absorbed into nmg_fuse.c.
+ * "meshing" here refers to the sorting of faceuses around an edge
+ * as two edges sharing the same end points (vertex structs) are fused.
+ *
+ */
+/** @} */
+
+#include "common.h"
+
+#include "nmg.h"
+#include "nmg_func.h"
+
+#include "bio.h"
+#include "bn/mat.h"
+#include "bn/plane.h"
+#include "bn/tol.h"
+#include "bu/log.h"
+
+#include <math.h>
+
+
+/**
+ * Determine if T lies within angle AB, such that A < T < B.
+ * The angle B is expected to be "more ccw" than A.
+ * Because of the wrap from 2pi to 0, B may have a smaller numeric value.
+ *
+ * Returns -
+ * -2 t is equal to a
+ * -1 t is equal to b
+ * 0 t is outside angle ab
+ * 1 t is inside angle ab
+ */
+HIDDEN int
+nmg_is_angle_in_wedge(double a, double b, double t)
+{
+    /* XXX What tolerance to use here (in radians)? */
+    if (NEAR_EQUAL(a, t, 1.0e-8)) return -2;
+    if (NEAR_EQUAL(b, t, 1.0e-8)) return -1;
+
+    /* If A==B, if T is not also equal, it's outside the wedge */
+    if (NEAR_EQUAL(a, b, 1.0e-8)) return 0;
+
+    if (b < a) {
+	/* B angle has wrapped past zero, add on 2pi */
+	if (t <= b) {
+	    /* Range is A..0, 0..B, and 0<t<B; so T is in wedge */
+	    return 1;
+	}
+	b += M_2PI;
+    }
+    if (NEAR_EQUAL(b, t, 1.0e-8)) return -1;
+
+    if (t < a) return 0;
+    if (t > b) return 0;
+    return 1;
+}
+
+
+/**
+ * Given two edgeuses with different edge geometry but
+ * running between the same two vertices,
+ * select the proper edge geometry to associate with.
+ *
+ * Really, there are 3 geometries to be compared here:
+ * the vector between the two endpoints of this edge,
+ * and the two edge_g structures.
+ * Rather than always taking eu2 or eu1,
+ * select the one that best fits this one edge.
+ *
+ * Consider fu1:
+ *		         B
+ *		         *
+ *		        /|
+ *		    eg2/ |
+ *		      /  |
+ *		    D/   |
+ *		    *    |
+ *		   /     |
+ *		A *-*----* C
+ *		    E eg1
+ *
+ * At the start of a face/face intersection, eg1 runs from A to C,
+ * and eg2 runs ADB.  The line of intersection with the other face
+ * (fu2, not drawn) lies along eg1.
+ * Assume that edge AC needs to be broken at E,
+ * where E is just a little more than tol->dist away from A.
+ * Existing point D is found because it *is* within tol->dist of E,
+ * thanks to the cosine of angle BAC.
+ * So, edge AC is broken on vertex D, and the intersection list
+ * contains vertexuses A, E, and C.
+ *
+ * Because D and E are the same point, fu1 has become a triangle with
+ * a little "spike" on the end.  If this is handled simply by re-homing
+ * edge AE to eg2, it may cause trouble, because eg1 now runs EC,
+ * but the geometry for eg1 runs AC.  If there are other vertices on
+ * edge eg1, the problem can not be resolved simply by recomputing the
+ * geometry of eg1.
+ * Since E (D) is within tolerance of eg1, it is not unreasonable
+ * just to leave eg1 alone.
+ *
+ * The issue boils down to selecting whether the existing eg1 or eg2
+ * best represents the direction of the little stub edge AD (shared with AE).
+ * In this case, eg2 is the correct choice, as AD (and AE) lie on line AB.
+ *
+ * It would be disastrous to force *all* of eg1 to use the edge geometry
+ * of eg2, as the two lines are very different.
+ */
+struct edge_g_lseg *
+nmg_pick_best_edge_g(struct edgeuse *eu1, struct edgeuse *eu2, const struct bn_tol *tol)
+{
+    NMG_CK_EDGEUSE(eu1);
+    NMG_CK_EDGEUSE(eu2);
+    BN_CK_TOL(tol);
+
+    NMG_CK_EDGE_G_LSEG(eu1->g.lseg_p);
+    NMG_CK_EDGE_G_LSEG(eu2->g.lseg_p);
+    if (eu2->g.lseg_p != eu1->g.lseg_p) {
+	vect_t dir;
+	vect_t dir_2;
+	vect_t dir_1;
+	fastf_t dot_2;
+	fastf_t dot_1;
+
+	VSUB2(dir, eu1->vu_p->v_p->vg_p->coord, eu1->eumate_p->vu_p->v_p->vg_p->coord);
+	VUNITIZE(dir);
+	VMOVE(dir_2, eu2->g.lseg_p->e_dir);
+	VUNITIZE(dir_2);
+	VMOVE(dir_1, eu1->g.lseg_p->e_dir);
+	VUNITIZE(dir_1);
+
+	dot_2 = fabs(VDOT(dir, dir_2));
+	dot_1 = fabs(VDOT(dir, dir_1));
+
+	/* Dot product of 1 means colinear.  Take largest dot. */
+	if (dot_2 > dot_1) {
+	    if (nmg_debug & DEBUG_BASIC) {
+		bu_log("nmg_pick_best_edge_g() Make eu1 use geometry of eu2, s.d=%g, d.d=%g\n",
+		       acos(dot_2)*RAD2DEG,
+		       acos(dot_1)*RAD2DEG);
+	    }
+	    return eu2->g.lseg_p;
+	} else {
+	    if (nmg_debug & DEBUG_BASIC) {
+		bu_log("nmg_pick_best_edge_g() Make eu2 use geometry of eu1, s.d=%g, d.d=%g\n",
+		       acos(dot_2)*RAD2DEG,
+		       acos(dot_1)*RAD2DEG);
+	    }
+	    return eu1->g.lseg_p;
+	}
+    }
+    return eu1->g.lseg_p;	/* both the same */
+}
+
+
+/**
+ * Make all the edgeuses around eu2's edge to refer to eu1's edge,
+ * taking care to organize them into the proper angular orientation,
+ * so that the attached faces are correctly arranged radially
+ * around the edge.
+ *
+ * This depends on both edges being part of face loops,
+ * with vertex and face geometry already associated.
+ *
+ * The two edgeuses being joined might well be from separate shells,
+ * so the issue of preserving (simple) faceuse orientation parity
+ * (SAME, OPPOSITE, OPPOSITE, SAME, ...)
+ * can't be used here -- that only applies to faceuses from the same shell.
+ *
+ * Some of the edgeuses around both edges may be wires.
+ *
+ * Call to nmg_check_radial at end has been deleted.
+ * Note that after two radial EU's have been joined
+ * a third cannot be joined to them without creating
+ * unclosed space that nmg_check_radial will find.
+ */
+void
+nmg_radial_join_eu(struct edgeuse *eu1, struct edgeuse *eu2, const struct bn_tol *tol)
+{
+
+    NMG_CK_EDGEUSE(eu1);
+    NMG_CK_EDGEUSE(eu1->radial_p);
+    NMG_CK_EDGEUSE(eu1->eumate_p);
+    NMG_CK_EDGEUSE(eu2);
+    NMG_CK_EDGEUSE(eu2->radial_p);
+    NMG_CK_EDGEUSE(eu2->eumate_p);
+    BN_CK_TOL(tol);
+
+    if (eu1->e_p == eu2->e_p) return;
+
+    if (!NMG_ARE_EUS_ADJACENT(eu1, eu2))
+	bu_bomb("nmg_radial_join_eu() edgeuses don't share vertices.\n");
+
+    if (eu1->vu_p->v_p == eu1->eumate_p->vu_p->v_p) bu_bomb("nmg_radial_join_eu(): 0 length edge (topology)\n");
+
+    if (bn_pt3_pt3_equal(eu1->vu_p->v_p->vg_p->coord,
+			 eu1->eumate_p->vu_p->v_p->vg_p->coord, tol))
+    {
+	bu_log("vertices should have been fused:\n");
+	bu_log("\tvertex %p (%.12f %.12f %.12f)\n",
+	       (void *)eu1->vu_p->v_p,
+	       V3ARGS(eu1->vu_p->v_p->vg_p->coord));
+	bu_log("\tvertex %p (%.12f %.12f %.12f)\n",
+	       (void *)eu1->eumate_p->vu_p->v_p,
+	       V3ARGS(eu1->eumate_p->vu_p->v_p->vg_p->coord));
+	bu_bomb("nmg_radial_join_eu(): 0 length edge (geometry)\n");
+    }
+
+#if 1
+    nmg_radial_join_eu_NEW(eu1, eu2, tol);
+    return;
+#else
+
+    /* Ensure faces are of same orientation, if both eu's have faces */
+    fu1 = nmg_find_fu_of_eu(eu1);
+    fu2 = nmg_find_fu_of_eu(eu2);
+    if (fu1 && fu2) {
+	if (fu1->orientation != fu2->orientation) {
+	    eu2 = eu2->eumate_p;
+	    fu2 = nmg_find_fu_of_eu(eu2);
+	    if (fu1->orientation != fu2->orientation)
+		bu_bomb("nmg_radial_join_eu(): Cannot find matching orientations for faceuses\n");
+	}
+    }
+
+    /* XXX This angle-based algorithm can't handle snurb faces! */
+    if (fu1 && fu1->f_p->g.magic_p && *fu1->f_p->g.magic_p == NMG_FACE_G_SNURB_MAGIC) return;
+    if (fu2 && fu2->f_p->g.magic_p && *fu2->f_p->g.magic_p == NMG_FACE_G_SNURB_MAGIC) return;
+
+    /* Construct local coordinate system for this edge,
+     * so all angles can be measured relative to a common reference.
+     */
+    nmg_eu_2vecs_perp(xvec, yvec, zvec, original_eu1, tol);
+
+    if (nmg_debug & DEBUG_MESH_EU) {
+	bu_log("nmg_radial_join_eu(eu1=%p, eu2=%p) e1=%p, e2=%p\n",
+	       (void *)eu1, (void *)eu2,
+	       (void *)eu1->e_p, (void *)eu2->e_p);
+	nmg_euprint("\tJoining", eu1);
+	nmg_euprint("\t     to", eu2);
+	bu_log("Faces around eu1:\n");
+	nmg_pr_fu_around_eu_vecs(eu1, xvec, yvec, zvec, tol);
+	bu_log("Faces around eu2:\n");
+	nmg_pr_fu_around_eu_vecs(eu2, xvec, yvec, zvec, tol);
+    }
+
+    best_eg = nmg_pick_best_edge_g(eu1, eu2, tol);
+
+    for (iteration1=0; eu2 && iteration1 < 10000; iteration1++) {
+	int code = 0;
+	struct edgeuse *first_eu1 = eu1;
+	int wire_skip = 0;
+	/* Resume where we left off from last eu2 insertion */
+
+	/* find a place to insert eu2 around eu1's edge */
+	for (iteration2=0; iteration2 < 10000; iteration2++) {
+	    struct faceuse *fur;
+
+	    abs1 = abs2 = absr = -M_2PI;
+
+	    eur = eu1->radial_p;
+	    NMG_CK_EDGEUSE(eur);
+
+	    fu2 = nmg_find_fu_of_eu(eu2);
+	    if (fu2 == (struct faceuse *)NULL) {
+		/* eu2 is a wire, it can go anywhere */
+		bu_log("eu2=%p is a wire, insert after eu1=%p\n", (void *)eu2, (void *)eu1);
+		goto insert;
+	    }
+	    fu1 = nmg_find_fu_of_eu(eu1);
+	    if (fu1 == (struct faceuse *)NULL) {
+		/* eu1 is a wire, skip on to real face eu */
+		bu_log("eu1=%p is a wire, skipping on\n", (void *)eu1);
+		wire_skip++;
+		goto cont;
+	    }
+	    fur = nmg_find_fu_of_eu(eur);
+	    while (fur == (struct faceuse *)NULL) {
+		/* eur is wire, advance eur */
+		bu_log("eur=%p is a wire, advancing to non-wire eur\n", (void *)eur);
+		eur = eur->eumate_p->radial_p;
+		wire_skip++;
+		if (eur == eu1->eumate_p) {
+		    bu_log("went all the way around\n");
+		    /* Went all the way around */
+		    goto insert;
+		}
+		fur = nmg_find_fu_of_eu(eur);
+	    }
+	    NMG_CK_FACEUSE(fu1);
+	    NMG_CK_FACEUSE(fu2);
+	    NMG_CK_FACEUSE(fur);
+
+	    /*
+	     * Can't just check for shared fg here,
+	     * the angle changes by +/- 180 degrees,
+	     * depending on which side of the eu the loop is on
+	     * along this edge.
+	     */
+	    abs1 = nmg_measure_fu_angle(eu1, xvec, yvec, zvec);
+	    abs2 = nmg_measure_fu_angle(eu2, xvec, yvec, zvec);
+	    absr = nmg_measure_fu_angle(eur, xvec, yvec, zvec);
+
+	    if (nmg_debug & DEBUG_MESH_EU) {
+		bu_log("  abs1=%g, abs2=%g, absr=%g\n",
+		       abs1*RAD2DEG,
+		       abs2*RAD2DEG,
+		       absr*RAD2DEG);
+	    }
+
+	    /* If abs1 == absr, warn about unfused faces, and skip. */
+	    if (NEAR_EQUAL(abs1, absr, 1.0e-8)) {
+		if (fu1->f_p->g.plane_p == fur->f_p->g.plane_p) {
+		    /* abs1 == absr, faces are fused, don't insert here. */
+		    if (nmg_debug & DEBUG_MESH_EU) {
+			bu_log("fu1 and fur share face geometry %p (flip1=%d, flip2=%d), skip\n",
+			       (void *)fu1->f_p->g.plane_p, fu1->f_p->flip, fur->f_p->flip);
+		    }
+		    goto cont;
+		}
+
+		bu_log("nmg_radial_join_eu: WARNING 2 faces should have been fused, may be ambiguous.\n  abs1=%e, absr=%e, asb2=%e\n",
+		       abs1*RAD2DEG, absr*RAD2DEG, abs2*RAD2DEG);
+		bu_log("  fu1=%p, f1=%p, f1->flip=%d, fg1=%p\n",
+		       (void *)fu1, (void *)fu1->f_p, fu1->f_p->flip, (void *)fu1->f_p->g.plane_p);
+		bu_log("  fu2=%p, f2=%p, f2->flip=%d, fg2=%p\n",
+		       (void *)fu2, (void *)fu2->f_p, fu2->f_p->flip, (void *)fu2->f_p->g.plane_p);
+		bu_log("  fur=%p, fr=%p, fr->flip=%d, fgr=%p\n",
+		       (void *)fur, (void *)fur->f_p, fur->f_p->flip, (void *)fur->f_p->g.plane_p);
+		PLPRINT("  fu1", fu1->f_p->g.plane_p->N);
+		PLPRINT("  fu2", fu2->f_p->g.plane_p->N);
+		PLPRINT("  fur", fur->f_p->g.plane_p->N);
+		{
+		    int debug = nmg_debug;
+
+		    if (nmg_two_face_fuse(fu1->f_p, fur->f_p, tol) == 0)
+			bu_bomb("faces didn't fuse?\n");
+		    nmg_debug = debug;
+		}
+		bu_log("  nmg_radial_join_eu() skipping this eu\n");
+		goto cont;
+	    }
+
+	    /*
+	     * If abs1 < abs2 < absr
+	     * (taking into account 360 wrap),
+	     * then insert face here.
+	     * Special handling if abs1==abs2 or abs2==absr.
+	     */
+	    code = nmg_is_angle_in_wedge(abs1, absr, abs2);
+	    if (nmg_debug & DEBUG_MESH_EU)
+		bu_log("    code=%d %s\n", code, (code!=0)?"INSERT_HERE":"skip");
+	    if (code > 0) break;
+	    if (code == -1) {
+		/* absr == abs2 */
+		break;
+	    }
+	    if (code <= -2) {
+		/* abs1 == abs2 */
+		break;
+	    }
+
+	cont:
+	    /* If eu1 is only one pair of edgeuses, done */
+	    if (eu1 == eur->eumate_p) break;
+	    eu1 = eur->eumate_p;
+	    if (eu1 == first_eu1) {
+		/* If all eu's were wires, here is fine */
+		if (wire_skip >= iteration2) break;
+		/* Nope, something bad happened */
+		bu_bomb("nmg_radial_join_eu():  went full circle, no face insertion point.\n");
+		break;
+	    }
+	}
+	if (iteration2 >= 10000) {
+	    bu_bomb("nmg_radial_join_eu: infinite loop (2)\n");
+	}
+
+	/* find the next use of the edge eu2 is on.  If eu2 and its
+	 * mate are the last uses of the edge, there will be no next
+	 * edgeuse to move. (Loop termination condition).
+	 */
+    insert:
+	nexteu = eu2->radial_p;
+	if (nexteu == eu2->eumate_p)
+	    nexteu = (struct edgeuse *)NULL;
+
+	/* because faces are always created with counter-clockwise
+	 * exterior loops and clockwise interior loops,
+	 * radial edgeuses IN THE SAME SHELL will never point in
+	 * the same direction or share the same vertex.  We thus make
+	 * sure that eu2 is an edgeuse which might be radial to eu1
+	 * XXX Need to look back for last eu IN THE SHELL OF eu2.
+	 * XXX Even this isn't good enough, as we may be inserting
+	 * XXX something new _after_ that last starting point.
+	 */
+	eus = eu1;
+	while (nmg_find_s_of_eu(eus) != nmg_find_s_of_eu(eu2)) {
+	    eus = eus->eumate_p->radial_p;
+	    if (eus == eu1) break;	/* full circle */
+	}
+
+	if (eu2->vu_p->v_p == eus->vu_p->v_p)
+	    eu2 = eu2->eumate_p;
+
+	if (nmg_debug & DEBUG_MESH_EU) {
+	    bu_log("  Inserting.  code=%d\n", code);
+	    bu_log("joining eu1=%p eu2=%p with abs1=%g, absr=%g\n",
+		   (void *)eu1, (void *)eu2,
+		   abs1*RAD2DEG, absr*RAD2DEG);
+	}
+
+	/*
+	 * Make eu2 radial to eu1.
+	 * This should insert eu2 between eu1 and eu1->radial_p
+	 * (which may be less far around than eur, but that's OK).
+	 * This does NOT change the edge geometry pointer.
+	 */
+	nmg_je(eu1, eu2);
+
+	if (nmg_debug & DEBUG_MESH_EU) {
+	    bu_log("After nmg_je(), faces around original_eu1 are:\n");
+	    nmg_pr_fu_around_eu_vecs(original_eu1, xvec, yvec, zvec, tol);
+	}
+
+	/* Proceed to the next source edgeuse */
+	eu2 = nexteu;
+    }
+    if (iteration1 >= 10000) bu_bomb("nmg_radial_join_eu:  infinite loop (1)\n");
+
+    NMG_CK_EDGEUSE(original_eu1);
+
+    /*
+     * Make another pass, ensuring that all edgeuses are using the
+     * "best_eg" line.
+     */
+    eu1 = original_eu1;
+    for (;;) {
+	if (eu1->g.lseg_p != best_eg) {
+	    nmg_use_edge_g(eu1, &best_eg->l.magic);
+	}
+
+	eu1 = eu1->eumate_p->radial_p;
+	if (eu1 == original_eu1) break;
+    }
+
+    if (nmg_debug & DEBUG_MESH_EU) bu_log("nmg_radial_join_eu: END\n");
+#endif
+}
+
+
+/**
+ * Actually do the work of meshing two faces.
+ * The two fu arguments may be the same, which causes the face to be
+ * meshed against itself.
+ *
+ * The return is the number of edges meshed.
+ */
+int
+nmg_mesh_two_faces(register struct faceuse *fu1, register struct faceuse *fu2, const struct bn_tol *tol)
+{
+    struct loopuse *lu1;
+    struct loopuse *lu2;
+    struct edgeuse *eu1;
+    struct edgeuse *eu2;
+    struct vertex *v1a, *v1b;
+    struct edge *e1;
+    pointp_t pt1, pt2;
+    int count = 0;
+
+    /* Visit all the loopuses in faceuse 1 */
+    for (BU_LIST_FOR(lu1, loopuse, &fu1->lu_hd)) {
+	/* Ignore self-loops */
+	if (BU_LIST_FIRST_MAGIC(&lu1->down_hd) != NMG_EDGEUSE_MAGIC)
+	    continue;
+
+	/* Visit all the edgeuses in loopuse1 */
+	for (BU_LIST_FOR(eu1, edgeuse, &lu1->down_hd)) {
+
+	    v1a = eu1->vu_p->v_p;
+	    v1b = eu1->eumate_p->vu_p->v_p;
+	    e1 = eu1->e_p;
+	    if (nmg_debug & DEBUG_MESH) {
+		pt1 = v1a->vg_p->coord;
+		pt2 = v1b->vg_p->coord;
+		bu_log("ref_e=%8p v:%8p--%8p (%g, %g, %g)->(%g, %g, %g)\n",
+		       (void *)e1, (void *)v1a, (void *)v1b,
+		       V3ARGS(pt1), V3ARGS(pt2));
+	    }
+
+	    /* Visit all the loopuses in faceuse2 */
+	    for (BU_LIST_FOR(lu2, loopuse, &fu2->lu_hd)) {
+		/* Ignore self-loops */
+		if (BU_LIST_FIRST_MAGIC(&lu2->down_hd) != NMG_EDGEUSE_MAGIC)
+		    continue;
+		/* Visit all the edgeuses in loopuse2 */
+		for (BU_LIST_FOR(eu2, edgeuse, &lu2->down_hd)) {
+		    if (nmg_debug & DEBUG_MESH) {
+			pt1 = eu2->vu_p->v_p->vg_p->coord;
+			pt2 = eu2->eumate_p->vu_p->v_p->vg_p->coord;
+			bu_log("\te:%8p v:%8p--%8p (%g, %g, %g)->(%g, %g, %g)\n",
+			       (void *)eu2->e_p,
+			       (void *)eu2->vu_p->v_p,
+			       (void *)eu2->eumate_p->vu_p->v_p,
+			       V3ARGS(pt1), V3ARGS(pt2));
+		    }
+
+		    /* See if already shared */
+		    if (eu2->e_p == e1) continue;
+		    if ((eu2->vu_p->v_p == v1a &&
+			 eu2->eumate_p->vu_p->v_p == v1b) ||
+			(eu2->eumate_p->vu_p->v_p == v1a &&
+			 eu2->vu_p->v_p == v1b)) {
+			nmg_radial_join_eu(eu1, eu2, tol);
+			count++;
+		    }
+		}
+	    }
+	}
+    }
+    return count;
+}
+
+
+/**
+ * Scan through all the edges of fu1 and fu2, ensuring that all
+ * edges involving the same vertex pair are indeed shared.
+ * This means worrying about merging ("meshing") all the faces in the
+ * proper radial orientation around the edge.
+ * XXX probably should return count;
+ */
+void
+nmg_mesh_faces(struct faceuse *fu1, struct faceuse *fu2, const struct bn_tol *tol)
+{
+    int count = 0;
+
+    NMG_CK_FACEUSE(fu1);
+    NMG_CK_FACEUSE(fu2);
+    BN_CK_TOL(tol);
+
+    if (nmg_debug & DEBUG_MESH_EU && nmg_debug & DEBUG_PLOTEM) {
+	nmg_pl_2fu("Before_mesh%d.plot3", fu1, fu2, 1);
+    }
+
+    if (nmg_debug & DEBUG_MESH_EU)
+	bu_log("meshing self (fu1 %8p)\n", (void *)fu1);
+    count += nmg_mesh_two_faces(fu1, fu1, tol);
+
+    if (nmg_debug & DEBUG_MESH_EU)
+	bu_log("meshing self (fu2 %8p)\n", (void *)fu2);
+    count += nmg_mesh_two_faces(fu2, fu2, tol);
+
+    if (nmg_debug & DEBUG_MESH_EU)
+	bu_log("meshing to other (fu1:%8p fu2:%8p)\n", (void *)fu1, (void *)fu2);
+    count += nmg_mesh_two_faces(fu1, fu2, tol);
+
+    if (nmg_debug & DEBUG_MESH_EU && nmg_debug & DEBUG_PLOTEM) {
+	nmg_pl_2fu("After_mesh%d.plot3", fu1, fu2, 1);
+    }
+}
+
+
+/**
+ * The return is the number of edges meshed.
+ */
+int
+nmg_mesh_face_shell(struct faceuse *fu1, struct shell *s, const struct bn_tol *tol)
+{
+    register struct faceuse *fu2;
+    int count = 0;
+
+    NMG_CK_FACEUSE(fu1);
+    NMG_CK_SHELL(s);
+    BN_CK_TOL(tol);
+
+    count += nmg_mesh_two_faces(fu1, fu1, tol);
+    for (BU_LIST_FOR(fu2, faceuse, &s->fu_hd)) {
+	NMG_CK_FACEUSE(fu2);
+	count += nmg_mesh_two_faces(fu2, fu2, tol);
+	count += nmg_mesh_two_faces(fu1, fu2, tol);
+    }
+    /* XXX What about wire edges in the shell? */
+    return count;
+}
+
+
+/**
+ * Mesh every edge in shell 1 with every edge in shell 2.
+ * The return is the number of edges meshed.
+ *
+ * Does not use nmg_mesh_face_shell() to keep face/self meshing
+ * to the absolute minimum necessary.
+ */
+int
+nmg_mesh_shell_shell(struct shell *s1, struct shell *s2, const struct bn_tol *tol)
+{
+    struct faceuse *fu1;
+    struct faceuse *fu2;
+    int count = 0;
+
+    NMG_CK_SHELL(s1);
+    NMG_CK_SHELL(s2);
+    BN_CK_TOL(tol);
+
+    nmg_region_v_unique(s1->r_p, tol);
+    nmg_region_v_unique(s2->r_p, tol);
+
+    /* First, mesh all faces of shell 2 with themselves */
+    for (BU_LIST_FOR(fu2, faceuse, &s2->fu_hd)) {
+	NMG_CK_FACEUSE(fu2);
+	count += nmg_mesh_two_faces(fu2, fu2, tol);
+    }
+
+    /* Visit every face in shell 1 */
+    for (BU_LIST_FOR(fu1, faceuse, &s1->fu_hd)) {
+	NMG_CK_FACEUSE(fu1);
+
+	/* First, mesh each face in shell 1 with itself */
+	count += nmg_mesh_two_faces(fu1, fu1, tol);
+
+	/* Visit every face in shell 2 */
+	for (BU_LIST_FOR(fu2, faceuse, &s2->fu_hd)) {
+	    NMG_CK_FACEUSE(fu2);
+	    count += nmg_mesh_two_faces(fu1, fu2, tol);
+	}
+    }
+
+    /* XXX What about wire edges in the shell? */
+
+    /* Visit every wire loop in shell 1 */
+
+    /* Visit every wire edge in shell 1 */
+
+    return count;
+}
+
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

src/libnmg/nmg_mirror.c

@@ -0,0 +1,172 @@
+/*                    N M G _ M I R R O R . C
+ * BRL-CAD
+ *
+ * Copyright (c) 2009-2016 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file primitives/nmg/nmg_mirror.c
+ *
+ * mirror support
+ *
+ */
+
+#include "common.h"
+
+#include "nmg.h"
+#include "nmg_func.h"
+
+#include "bn/mat.h"
+#include "bn/tol.h"
+#include "bu/log.h"
+
+
+/**
+ * Given a pointer to an internal GED database object, mirror the
+ * object's values about the given transformation matrix.
+ */
+HIDDEN int
+rt_nmg_mirror(struct model *nmg, register const plane_t plane)
+{
+    mat_t mirmat;
+    mat_t rmat;
+    mat_t temp;
+    vect_t nvec;
+    vect_t xvec;
+    vect_t mirror_dir;
+    point_t mirror_pt;
+    fastf_t ang;
+
+    int i;
+    struct nmgregion *r;
+    struct shell *s;
+    struct bu_ptbl table;
+    struct vertex *v;
+
+    static point_t origin = {0.0, 0.0, 0.0};
+
+    static const struct bn_tol tol = {
+	BN_TOL_MAGIC, 0.0005, 0.0005*0.0005, 1e-6, 1-1e-6
+    };
+
+    NMG_CK_MODEL(nmg);
+
+    MAT_IDN(mirmat);
+
+    VMOVE(mirror_dir, plane);
+    VSCALE(mirror_pt, plane, plane[W]);
+
+    /* Build mirror transform matrix, for those who need it. */
+    /* First, perform a mirror down the X axis */
+    mirmat[0] = -1.0;
+
+    /* Create the rotation matrix */
+    VSET(xvec, 1, 0, 0);
+    VCROSS(nvec, xvec, mirror_dir);
+    VUNITIZE(nvec);
+    ang = -acos(VDOT(xvec, mirror_dir));
+    bn_mat_arb_rot(rmat, origin, nvec, ang*2.0);
+
+    /* Add the rotation to mirmat */
+    MAT_COPY(temp, mirmat);
+    bn_mat_mul(mirmat, temp, rmat);
+
+    /* Add the translation to mirmat */
+    mirmat[3 + X*4] += mirror_pt[X] * mirror_dir[X];
+    mirmat[3 + Y*4] += mirror_pt[Y] * mirror_dir[Y];
+    mirmat[3 + Z*4] += mirror_pt[Z] * mirror_dir[Z];
+
+    /* move every vertex */
+    nmg_vertex_tabulate(&table, &nmg->magic);
+    for (i=0; i<BU_PTBL_END(&table); i++) {
+	point_t pt;
+
+	v = (struct vertex *)BU_PTBL_GET(&table, i);
+	NMG_CK_VERTEX(v);
+
+	VMOVE(pt, v->vg_p->coord);
+	MAT4X3PNT(v->vg_p->coord, mirmat, pt);
+    }
+
+    bu_ptbl_reset(&table);
+
+    nmg_face_tabulate(&table, &nmg->magic);
+    for (i=0; i<BU_PTBL_END(&table); i++) {
+	struct face *f;
+
+	f = (struct face *)BU_PTBL_GET(&table, i);
+	NMG_CK_FACE(f);
+
+	if (!f->g.magic_p)
+	    continue;
+
+	if (*f->g.magic_p != NMG_FACE_G_PLANE_MAGIC) {
+	    bu_log("Sorry, can only mirror NMG solids with planar faces\n");
+	    bu_ptbl_free(&table);
+	    return 1;
+	}
+    }
+
+    for (BU_LIST_FOR (r, nmgregion, &nmg->r_hd)) {
+	for (BU_LIST_FOR (s, shell, &r->s_hd)) {
+	    nmg_invert_shell(s);
+	}
+    }
+
+
+    for (i=0; i<BU_PTBL_END(&table); i++) {
+	struct face *f;
+	struct faceuse *fu;
+
+	f = (struct face *)BU_PTBL_GET(&table, i);
+	NMG_CK_FACE(f);
+
+	fu = f->fu_p;
+	if (fu->orientation != OT_SAME) {
+	    fu = fu->fumate_p;
+	}
+	if (fu->orientation != OT_SAME) {
+	    bu_log("ERROR: Unexpected NMG face orientation\n");
+	    bu_ptbl_free(&table);
+	    return 2;
+	}
+
+	if (nmg_calc_face_g(fu)) {
+	    bu_log("ERROR: Unable to calculate NMG faces for mirroring\n");
+	    bu_ptbl_free(&table);
+	    return 3;
+	}
+    }
+
+    bu_ptbl_free(&table);
+
+    /* FIXME: why do we need to rebound?  should mirroring really
+     * require a tolerance??
+     */
+    nmg_rebound(nmg, &tol);
+
+    return 0;
+}
+
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

src/libnmg/nmg_visit.c

@@ -0,0 +1,304 @@
+/*                     N M G _ V I S I T . C
+ * BRL-CAD
+ *
+ * Copyright (c) 1993-2016 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @addtogroup nmg */
+/** @{ */
+/** @file primitives/nmg/nmg_visit.c
+ *
+ * A generalized, object-oriented subroutine family to
+ * visit all the data structures "below" a given structure.
+ *
+ * The caller provides a pointer to the structure to start at,
+ * a table of "handlers" for each kind of structure,
+ * and a generic pointer for private state which will be sent along
+ * to the user's handlers.
+ * For non-leaf structures, there are two handlers, one called
+ * before any recursion starts, and the other called when
+ * recursion is finished.  Either or both may be omitted.
+ *
+ */
+/** @} */
+
+#include "common.h"
+
+#include "nmg.h"
+
+#include "bio.h"
+#include "bu/log.h"
+
+#include <math.h>
+
+
+void
+nmg_visit_vertex(struct vertex *v, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+    NMG_CK_VERTEX(v);
+
+    if (htab->vis_vertex) htab->vis_vertex((uint32_t *)v, state, 0);
+
+    if (htab->vis_vertex_g && v->vg_p)
+	htab->vis_vertex_g((uint32_t *)v->vg_p, state, 0);
+}
+
+
+void
+nmg_visit_vertexuse(struct vertexuse *vu, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+    NMG_CK_VERTEXUSE(vu);
+
+    if (htab->bef_vertexuse) htab->bef_vertexuse((uint32_t *)vu, state, 0);
+
+    nmg_visit_vertex(vu->v_p, htab, state);
+
+    if (htab->vis_vertexuse_a && vu->a.magic_p)
+	htab->vis_vertexuse_a((uint32_t *)vu->a.magic_p, state, 0);
+
+    if (htab->aft_vertexuse) htab->aft_vertexuse((uint32_t *)vu, state, 1);
+}
+
+
+void
+nmg_visit_edge(struct edge *e, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+    NMG_CK_EDGE(e);
+
+    if (htab->vis_edge) htab->vis_edge((uint32_t *)e, state, 0);
+}
+
+
+void
+nmg_visit_edgeuse(struct edgeuse *eu, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+    NMG_CK_EDGEUSE(eu);
+
+    if (htab->bef_edgeuse) htab->bef_edgeuse((uint32_t *)eu, state, 0);
+
+    nmg_visit_vertexuse(eu->vu_p, htab, state);
+    nmg_visit_edge(eu->e_p, htab, state);
+
+    if (htab->vis_edge_g && eu->g.magic_p)
+	htab->vis_edge_g((uint32_t *)eu->g.magic_p, state, 0);
+
+    if (htab->aft_edgeuse) htab->aft_edgeuse((uint32_t *)eu, state, 1);
+}
+
+
+void
+nmg_visit_loop(struct loop *l, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+    NMG_CK_LOOP(l);
+
+    if (htab->vis_loop) htab->vis_loop((uint32_t *)l, state, 0);
+
+    if (htab->vis_loop_g && l->lg_p)
+	htab->vis_loop_g((uint32_t *)l->lg_p, state, 0);
+}
+
+
+void
+nmg_visit_loopuse(struct loopuse *lu, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+    NMG_CK_LOOPUSE(lu);
+
+    if (htab->bef_loopuse) htab->bef_loopuse((uint32_t *)lu, state, 0);
+
+    if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
+	struct vertexuse *vu;
+	vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
+	nmg_visit_vertexuse(vu, htab, state);
+    } else {
+	struct edgeuse *eu;
+	for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
+	    nmg_visit_edgeuse(eu, htab, state);
+	}
+    }
+    nmg_visit_loop(lu->l_p, htab, state);
+
+    if (htab->aft_loopuse) htab->aft_loopuse((uint32_t *)lu, state, 1);
+}
+
+
+void
+nmg_visit_face(struct face *f, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+
+    if (htab->vis_face) htab->vis_face((uint32_t *)f, state, 0);
+
+    if (htab->vis_face_g && f->g.plane_p)
+	htab->vis_face_g((uint32_t *)f->g.plane_p, state, 0);
+}
+
+
+void
+nmg_visit_faceuse(struct faceuse *fu, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+    struct loopuse *lu;
+
+    NMG_CK_FACEUSE(fu);
+
+    if (htab->bef_faceuse) htab->bef_faceuse((uint32_t *)fu, state, 0);
+
+    for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
+	nmg_visit_loopuse(lu, htab, state);
+    }
+
+    nmg_visit_face(fu->f_p, htab, state);
+
+    if (htab->aft_faceuse) htab->aft_faceuse((uint32_t *)fu, state, 1);
+}
+
+
+void
+nmg_visit_shell(struct shell *s, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+    struct faceuse *fu;
+    struct loopuse *lu;
+    struct edgeuse *eu;
+
+    NMG_CK_SHELL(s);
+
+    if (htab->bef_shell) htab->bef_shell((uint32_t *)s, state, 0);
+
+    for (BU_LIST_FOR(fu, faceuse, &s->fu_hd)) {
+	nmg_visit_faceuse(fu, htab, state);
+    }
+    for (BU_LIST_FOR(lu, loopuse, &s->lu_hd)) {
+	nmg_visit_loopuse(lu, htab, state);
+    }
+    for (BU_LIST_FOR(eu, edgeuse, &s->eu_hd)) {
+	nmg_visit_edgeuse(eu, htab, state);
+    }
+    if (s->vu_p) nmg_visit_vertexuse(s->vu_p, htab, state);
+    if (htab->vis_shell_a && s->sa_p)
+	htab->vis_shell_a((uint32_t *)s->sa_p, state, 0);
+
+    if (htab->aft_shell) htab->aft_shell((uint32_t *)s, state, 1);
+}
+
+
+void
+nmg_visit_region(struct nmgregion *r, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+    struct shell *s;
+
+    NMG_CK_REGION(r);
+
+    if (htab->bef_region) htab->bef_region((uint32_t *)r, state, 0);
+
+    for (BU_LIST_FOR(s, shell, &r->s_hd)) {
+	nmg_visit_shell(s, htab, state);
+    }
+    if (htab->vis_region_a && r->ra_p)
+	htab->vis_region_a((uint32_t *)r->ra_p, state, 0);
+
+    if (htab->aft_region) htab->aft_region((uint32_t *)r, state, 1);
+}
+void
+nmg_visit_model(struct model *model, const struct nmg_visit_handlers *htab, void *state)
+
+
+/* Handler's private state */
+{
+    struct nmgregion *r;
+
+    NMG_CK_MODEL(model);
+
+    if (htab->bef_model) htab->bef_model((uint32_t *)model, state, 0);
+
+    for (BU_LIST_FOR(r, nmgregion, &model->r_hd)) {
+	nmg_visit_region(r, htab, state);
+    }
+
+    if (htab->aft_model) htab->aft_model((uint32_t *)model, state, 1);
+}
+
+
+void
+nmg_visit(const uint32_t *magicp, const struct nmg_visit_handlers *htab, void *state)
+/* Handler's private state */
+{
+    switch (*magicp) {
+	default:
+	    bu_log("nmg_visit() Can't visit %s directly\n", bu_identify_magic(*magicp));
+	    bu_bomb("nmg_visit()\n");
+	    /* NOTREACHED */
+	case NMG_MODEL_MAGIC:
+	    nmg_visit_model((struct model *)magicp, htab, state);
+	    break;
+	case NMG_REGION_MAGIC:
+	    nmg_visit_region((struct nmgregion *)magicp, htab, state);
+	    break;
+	case NMG_SHELL_MAGIC:
+	    nmg_visit_shell((struct shell *)magicp, htab, state);
+	    break;
+	case NMG_FACEUSE_MAGIC:
+	    nmg_visit_faceuse((struct faceuse *)magicp, htab, state);
+	    break;
+	case NMG_LOOPUSE_MAGIC:
+	    nmg_visit_loopuse((struct loopuse *)magicp, htab, state);
+	    break;
+	case NMG_EDGEUSE_MAGIC:
+	    nmg_visit_edgeuse((struct edgeuse *)magicp, htab, state);
+	    break;
+	case NMG_VERTEXUSE_MAGIC:
+	    nmg_visit_vertexuse((struct vertexuse *)magicp, htab, state);
+	    break;
+    }
+}
+
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */