geometry: check overlap on final segments and class junction nodes

src/c_geometry.h

@@ -27,9 +27,22 @@
 
 class nec_context;
 #include "nec_output.h"
-#include "nec_wire.h"
 
 
+/*! \brief One segment center found inside the volume of a segment that it does
+ *  not meet at a shared node.
+ *
+ *  The distance and the radius are the measurement that produced the finding,
+ *  so a caller can judge how far inside the conductor the center falls.
+ */
+struct nec_overlap_finding {
+  int64_t   inside_index;
+  int64_t   container_index;
+  int       inside_tag;
+  int       container_tag;
+  nec_float distance;
+  nec_float container_radius;
+};
 
 /*! \brief A Class describing the antenna geometry
  * \file c_geometry.h
@@ -134,6 +147,18 @@ class c_geometry
   	 */
   	void set_intersection_check(bool enable) { _check_intersections = enable; }
 
+  	/*! \brief Choose what an overlap finding does.
+  	 *  \param fatal true to reject the geometry on the first finding (the
+  	 *  default), false to record every finding and let the solve proceed.
+  	 */
+  	void set_intersection_fatal(bool fatal) { _overlap_is_fatal = fatal; }
+
+  	/*! \brief The overlap findings the last geometry check produced.
+  	 *  \return One record per segment center found inside the volume of a
+  	 *  segment it does not meet at a node.
+  	 */
+  	const std::vector<nec_overlap_finding>& overlap_findings() const { return m_overlap_findings; }
+
   	/*! \brief Geometry is complete
   	 *  \exception nec_exception* If there is an error with the geometry.
   	 */
@@ -208,6 +233,17 @@ class c_geometry
   void divide_patch( int nx );
 
   void connect_segments( int ignd );
+
+  /*! \brief Record every segment center that lies inside the volume of a
+   *  segment it is not connected to, then apply the overlap policy to the
+   *  findings.
+   *  \exception nec_exception* Under a fatal policy, if two unconnected
+   *  segments overlap.
+   */
+  void check_segment_intersections();
+
+  void collect_overlap_findings();
+
   void build_connections( int ignd );
   void resolve_junctions();
 
@@ -232,8 +268,8 @@ class c_geometry
   	nec_context* m_context;
   	nec_output_file* m_output;
   	bool _check_intersections;
-
-  	std::vector<nec_wire> m_wires;
+  	bool _overlap_is_fatal;
+  	std::vector<nec_overlap_finding> m_overlap_findings;
 
   void reflect_plane(int sym_plane, int& tag_increment);
 

src/misc.cpp

@@ -24,6 +24,7 @@ void usage(void)
       "\n       -s: print result summary to standard output."
       "\n       -c: print results in comma-separated-value (CSV) format,"
       "\n           this options is used in conjunction with (-s) above."
+      "\n       -w: report a geometry overlap as a warning and continue."
       "\n       -h: print this usage information and exit."
       "\n       -v: print nec2++ version number and exit.\n";
 

src/nec2cpp.cpp

@@ -180,7 +180,7 @@ int nec_main( int argc, char **argv, nec_output_file& s_output )
 	nec_context s_context;
 
 	/* process command line options */
-	while( (option = XGetopt(argc, argv, "i:o:hvscxgb") ) != -1 )
+	while( (option = XGetopt(argc, argv, "i:o:hvscxgbw") ) != -1 )
 	{
 		switch( option )
 		{
@@ -208,6 +208,10 @@ int nec_main( int argc, char **argv, nec_output_file& s_output )
 			s_context.set_results_format(RESULT_FORMAT_XML);
 			break;
 
+		case 'w': /* report a geometry overlap as a warning and continue */
+			s_context.get_geometry()->set_intersection_fatal(false);
+			break;
+
 		case 'h' : /* print usage and exit */
 			usage();
 			exit(0);

src/nec_context.cpp

@@ -96,6 +96,8 @@ void nec_context::initialize() {
   iptaqf=0;
   iptaqt=0;
 
+  step_warning_issued = false;
+
   init_voltage_sources();
   m_geometry->set_context(this);
 }
@@ -5832,7 +5834,6 @@ void nec_context::rom2( nec_float a, nec_float b, complex_array& sum, nec_float
   ASSERT(sum.size() == 9);
 
   bool recalculate_fields = true;
-  static bool step_warning_issued = false;
 
   int nts = 4, nx = 1, n = 9;
 

src/nec_context.h

@@ -922,6 +922,11 @@ class nec_context
 
 private:
 
+  /*! \brief Set once the ROM2 step-size explanation has been emitted, so this
+   *  context explains it one time rather than once per limited integration.
+   */
+  bool step_warning_issued;
+
   /*! \brief A private convenience function called by ne_card() and nh_card()
   */
   void ne_nh_card(int in_nfeh, int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6);

src/nec_context_tb.cpp

@@ -211,6 +211,138 @@ TEST_CASE( "Optional intersection check bypass", "[intersection_check]") {
     REQUIRE_NOTHROW(nec.geometry_complete(0));
 }
 
+TEST_CASE( "Crossed wires sharing a node are accepted", "[segment_intersection]") {
+    // Turnstile geometry: two dipoles crossing at the origin, each contributing
+    // a segment end to the shared node. NEC-2 places no restriction on the angle
+    // between connected wires, and every segment center clears the crossing
+    // wire's volume by several radii.
+    nec_context nec;
+    nec.initialize();
+
+    c_geometry* geo = nec.get_geometry();
+    geo->wire(1, 24, 0.0, 0.517, 0.0, 0.0, -0.517, 0.0, 0.006, 1.0, 1.0);
+    geo->wire(2, 24, -0.517, 0.0, 0.0, 0.517, 0.0, 0.0, 0.006, 1.0, 1.0);
+    REQUIRE_NOTHROW(nec.geometry_complete(0));
+}
+
+TEST_CASE( "Translated wire is tested at its final position", "[segment_intersection]") {
+    // Tower geometry: a leg is built from the ground up, translated by its own
+    // base height, and a second leg then fills the span it vacated. The two legs
+    // meet end to end and share no volume in the final segment arrays.
+    nec_context nec;
+    nec.initialize();
+
+    c_geometry* geo = nec.get_geometry();
+    geo->wire(3, 20, 1.249, 0.0, 0.0, 1.249, 0.0, 11.24, 0.05, 1.0, 1.0);
+    geo->move(0.0, 0.0, 0.0, 0.0, 0.0, 1.249, 0, 0, 0);
+    geo->wire(1, 3, 1.249, 0.0, 0.0, 1.249, 0.0, 1.249, 0.05, 1.0, 1.0);
+    REQUIRE_NOTHROW(nec.geometry_complete(0));
+}
+
+TEST_CASE( "Tapered wire is tested at its per-segment radius", "[segment_intersection]") {
+    // Biconical geometry: a tapered cone meeting a short centre segment. The
+    // cone's widest radius belongs to its far segment, and the segment actually
+    // adjoining the centre is thin enough to clear it.
+    nec_context nec;
+    nec.initialize();
+
+    c_geometry* geo = nec.get_geometry();
+    geo->wire(1, 5, 0.0, 0.0, -1.7831, 0.0, 0.0, -0.1024, 0.1514, 1.0, 0.6399);
+    geo->wire(2, 1, 0.0, 0.0, -0.1024, 0.0, 0.0, 0.1024, 0.0254, 1.0, 1.0);
+    REQUIRE_NOTHROW(nec.geometry_complete(0));
+}
+
+TEST_CASE( "Coincident wires are rejected", "[segment_intersection]") {
+    // Two wires occupying the same space: every segment center of one lies on
+    // the axis of a segment of the other.
+    nec_context nec;
+    nec.initialize();
+
+    c_geometry* geo = nec.get_geometry();
+    geo->wire(1, 5, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.001, 1.0, 1.0);
+    geo->wire(2, 5, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.001, 1.0, 1.0);
+    REQUIRE_THROWS(nec.geometry_complete(0));
+}
+
+TEST_CASE( "Wire penetrating another away from any node is rejected", "[segment_intersection]") {
+    // Two single-segment wires crossing at their centres, sharing no end, so the
+    // centre of each falls inside the volume of the other.
+    nec_context nec;
+    nec.initialize();
+
+    c_geometry* geo = nec.get_geometry();
+    geo->wire(1, 1, -0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.01, 1.0, 1.0);
+    geo->wire(2, 1, 0.0, -0.5, 0.0, 0.0, 0.5, 0.0, 0.01, 1.0, 1.0);
+    REQUIRE_THROWS(nec.geometry_complete(0));
+}
+
+TEST_CASE( "Six wires meeting at one node are accepted", "[segment_junction]") {
+    // Ground plane vertical on a mast: four radials, a radiator, and the mast
+    // all contribute an end to the node at the origin. The connection record
+    // links those six ends as a cycle, so the radials are not adjacent to the
+    // mast in it, while every one of them is joined to it in the geometry.
+    nec_context nec;
+    nec.initialize();
+
+    c_geometry* geo = nec.get_geometry();
+    geo->wire(1, 13, 0.0, 0.0, 0.0, -0.34, 0.0, -0.34, 0.0075, 1.0, 1.0);
+    geo->wire(1, 13, 0.0, 0.0, 0.0, 0.34, 0.0, -0.34, 0.0075, 1.0, 1.0);
+    geo->wire(1, 13, 0.0, 0.0, 0.0, 0.0, -0.34, -0.34, 0.0075, 1.0, 1.0);
+    geo->wire(1, 13, 0.0, 0.0, 0.0, 0.0, 0.34, -0.34, 0.0075, 1.0, 1.0);
+    geo->wire(2, 13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0075, 1.0, 1.0);
+    geo->wire(3, 75, 0.0, 0.0, 0.0, 0.0, 0.0, -3.0, 0.025, 1.0, 1.0);
+    REQUIRE_NOTHROW(nec.geometry_complete(0));
+    REQUIRE(geo->overlap_findings().empty());
+}
+
+TEST_CASE( "Legs meeting on the ground plane are accepted", "[segment_junction]") {
+    // Two legs leaving one point on the ground plane 45 degrees apart, each
+    // first center inside the other's volume and every later center clear of
+    // it. build_connections() records a ground contact rather than a neighbour
+    // at that end, so the junction is known from the coordinates alone.
+    nec_context nec;
+    nec.initialize();
+
+    c_geometry* geo = nec.get_geometry();
+    geo->wire(1, 10, 0.0, 0.0, 0.0, 0.15307, 0.0, 0.36955, 0.02, 1.0, 1.0);
+    geo->wire(2, 10, 0.0, 0.0, 0.0, -0.15307, 0.0, 0.36955, 0.02, 1.0, 1.0);
+    REQUIRE_NOTHROW(nec.geometry_complete(1));
+    REQUIRE(geo->overlap_findings().empty());
+}
+
+TEST_CASE( "Wire intruding past a shared node is rejected", "[segment_junction]") {
+    // A feed leaving the hub of a radial at 20 degrees. Its first center clears
+    // the radial segment it shares the hub with and comes to rest inside the
+    // next one along, which it meets at no node.
+    nec_context nec;
+    nec.initialize();
+
+    c_geometry* geo = nec.get_geometry();
+    geo->wire(2, 10, 0.0, 0.0, 0.0, 0.045, 0.0, 0.0, 0.003, 1.0, 1.0);
+    geo->wire(3, 5, 0.0, 0.0, 0.0, 0.07518, 0.0, 0.02736, 0.001, 1.0, 1.0);
+    REQUIRE_THROWS(nec.geometry_complete(0));
+}
+
+TEST_CASE( "Warning policy records the overlap and continues", "[segment_junction]") {
+    // The geometry of the coincident-wire rejection, admitted under a policy
+    // that reports rather than rejects. Every center of one wire lies on the
+    // axis of the other, so both directions of every pair are recorded.
+    nec_context nec;
+    nec.initialize();
+
+    c_geometry* geo = nec.get_geometry();
+    geo->set_intersection_fatal(false);
+    geo->wire(1, 5, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.001, 1.0, 1.0);
+    geo->wire(2, 5, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.001, 1.0, 1.0);
+    REQUIRE_NOTHROW(nec.geometry_complete(0));
+
+    REQUIRE(false == geo->overlap_findings().empty());
+    const nec_overlap_finding& first = geo->overlap_findings().front();
+    REQUIRE(first.inside_tag == 1);
+    REQUIRE(first.container_tag == 2);
+    REQUIRE(first.distance <= first.container_radius);
+}
+
 TEST_CASE( "Helix rejects invalid segment_count", "[helix]") {
     // Regression test for #48: helix with segment_count < 1 must throw
     // rather than silently returning with no wires.

src/nec_wire.h

@@ -94,6 +94,16 @@ std::vector<nec_wire> intersect(nec_wire& b)
 \return true if the point x is inside the wire
 */
 	bool intersect(nec_3vector& b0)
+	{
+		return axis_distance(b0) <= radius;
+	}
+
+/*!\brief Measure how far a point lies from the wire axis
+\param b0 The point being measured.
+\return The distance from b0 to the nearest point of the axis, which is clamped
+	to the wire's own extent so a point beyond an end measures from that end.
+*/
+	nec_float axis_distance(const nec_3vector& b0) const
 	{
 		nec_float a0x = x0(0); nec_float a0y = x0(1); nec_float a0z = x0(2);
 		nec_float a1x = x1(0); nec_float a1y = x1(1); nec_float a1z = x1(2);
@@ -159,9 +169,8 @@ print solution
 
   nec_3vector a_pt = parametrize(sa);
 
-	if (distance(a_pt, b0) > radius) return false;
-	return true;
-
+	return distance(a_pt, b0);
+
 }
 
 	bool similar(nec_wire& b)