From 94774de1d385c719e8217743a9ce34216175842b Mon Sep 17 00:00:00 2001 From: Eric Wheeler Date: Wed, 5 Aug 2026 18:53:23 -0700 Subject: [PATCH 1/3] geometry: derive intersection check from final segment arrays The pre-solve wire intersection check read m_wires, a shadow list of the original GW card inputs, instead of the segment arrays the solver actually uses. That shadow list carried no segmentation topology, was never updated by GM, GX, GR, or GS, and stored one radius per card even when GC tapered each segment along its length. Valid NEC-2 models were rejected at crossed nodes, translated wires, and tapered cones, while some genuinely overlapping models went undetected because parallel axes drove the whole-wire solve into a sentinel return. - delete m_wires and its writes in wire(), so geometry has one owner: the segment arrays populated after connect_segments() establishes topology - add check_segment_intersections(), sweeping final segment centers, lengths, radii, and direction cosines along the structure's widest axis, bounded by each segment's own half-length-plus-radius reach - exclude segment pairs NEC already recorded as joined via icon1/icon2, so a shared node reads as a connection rather than an intersection - call the new check from geometry_complete() after the SEGMENT DATA ERROR gate, so no zero-length axis reaches the point-in-cylinder solve - add permanent regressions for crossed wires sharing a node, a wire translated by move(), a tapered chain, coincident wires, and interior penetration away from any node Signed-off-by: Eric Wheeler --- src/c_geometry.cpp | 324 +++++++++++++++++++++++++++-------------- src/c_geometry.h | 10 +- src/nec_context_tb.cpp | 65 +++++++++ src/nec_wire.h | 15 +- 4 files changed, 299 insertions(+), 115 deletions(-) diff --git a/src/c_geometry.cpp b/src/c_geometry.cpp index 9f3899bf..6e6d216c 100644 --- a/src/c_geometry.cpp +++ b/src/c_geometry.cpp @@ -20,11 +20,25 @@ #include "nec_context.h" #include "nec_exception.h" +#include "nec_wire.h" #include #include #include +/*!\brief One final segment as the intersection sweep reads it. + + Carries the segment's own index so the sorted sweep still names it, its final + geometry as a wire, its center, and the segments its two ends are joined to. +*/ +struct segment_view { + int64_t index; + nec_wire body; + nec_3vector midpoint; + int64_t linked_start; + int64_t linked_end; +}; + /** * geometry_field_separator - Determine whether a character separates fields * @character: Character read from a geometry card @@ -445,42 +459,14 @@ void c_geometry::parse_geometry_error(const geometry_parse_state& st) throw nec_exception("GEOMETRY DATA CARD ERROR"); } -#include "nec_wire.h" /** - We have finished with the geometry description, now connect + We have finished with the geometry description, now connect things up. */ void c_geometry::geometry_complete(nec_context* in_context, int gpflag) { if (0 == np + mp) throw nec_exception("Geometry has no wires or patches."); - - /* Check to see whether any wires intersect with one another */ -if (_check_intersections) -{ - for (size_t i=0; i j) - { - nec_wire b = m_wires[j]; - vector wires = a.intersect(b); - if (wires.size() > 2) - { - nec_exception nex("GEOMETRY DATA ERROR -- WIRE #"); - nex.append(j+1); - nex.append(" (TAG ID #"); nex.append(b.tag_id()); - nex.append(") INTERSECTS WIRE #"); - nex.append(i+1); - nex.append(" (TAG ID #"); nex.append(a.tag_id()); nex.append(")"); - throw nex; - } - } - } - } -} // now proceed and complete the geometry setup... // Check here that patches form a closed surfaceAntennaInput @@ -566,6 +552,10 @@ if (_check_intersections) } /* for( i = 0; i < n_segments; i++ ) */ } /* if ( n_segments != 0) */ + // Segment centers, lengths, direction cosines, and junction records are all + // final here, so the overlap rule is applied to the geometry the solver uses. + check_segment_intersections(); + if ( m != 0) { m_output->nec_printf( "\n\n\n" @@ -677,67 +667,8 @@ void c_geometry::wire( int tag_id, int segment_count, nec_float xw1, nec_float y rd=1.0; } - /* - There is no restriction on the angle between two wires, but accuracy will be lost if the center of a segment falls within the - volume of the wire the segment connects to. The risk of this reduces as the angle between wires approaches 180 degrees. - - Wires which intersect away from their ends are not connected, but errors will occur if one wire occupies the space of another one. For accuracy, separate wire centers by several radii of the largest wire. - */ - // check that none of the existing wires intersect with the midpoint - // of the first and last segment - nec_3vector wire_start(xw1,yw1,zw1); - nec_3vector wire_end(xw2,yw2,zw2); - - if (_check_intersections) - { - nec_3vector seg_midpoint(wire_start + (dx/2)*delz); - nec_3vector end_seg_midpoint = wire_end - (dx*delz / 2); - /* Check to see whether any wires intersect with the segment_midpoint */ - for (uint32_t i=0; i(connection) + : static_cast(connection); + bool joins_segment = (0 != connection) && (magnitude <= PCHCON); + return joins_segment ? (magnitude - 1) : -1; +} + +/*! \brief Report whether NEC recorded these two segments as sharing an end. + \param a One segment view. + \param b The other segment view. + \return true if either segment records the other at either of its ends. + + build_connections() stops at the first contact found for each end, so a + junction of three or more segments is recorded as a chain rather than a + clique. Reading both ends of both segments covers every pair that chain + records. +*/ +static bool segments_joined(const segment_view& a, const segment_view& b) +{ + return (a.linked_start == b.index) || (a.linked_end == b.index) || + (b.linked_start == a.index) || (b.linked_end == a.index); +} + +/*! \brief Report one segment center inside another segment. + \param inside The segment whose center falls within the other's volume. + \param container The segment whose volume contains that center. + \exception nec_exception* Always thrown, naming both segments. +*/ +static void throw_segment_overlap(const segment_view& inside, + const segment_view& container) +{ + nec_exception nex("GEOMETRY DATA ERROR -- SEGMENT #"); + nex.append(inside.index + 1); + nex.append(" (TAG ID #"); nex.append(inside.body.tag_id()); + nex.append(") MIDPOINT LIES WITHIN SEGMENT #"); + nex.append(container.index + 1); + nex.append(" (TAG ID #"); nex.append(container.body.tag_id()); + nex.append(")"); + throw nex; +} + +/*! \brief Find the axis over which the structure spreads furthest. + \param views One view per final segment. + \return The coordinate index to sweep along. + + Sweeping along this axis leaves the ordered window admitting the fewest + candidate pairs. +*/ +static Eigen::Index widest_extent_axis(const std::vector& views) +{ + nec_3vector low = views[0].midpoint; + nec_3vector high = low; + for (size_t i = 1; i < views.size(); i++) + { + low = low.cwiseMin(views[i].midpoint); + high = high.cwiseMax(views[i].midpoint); + } + + Eigen::Index axis = 0; + (high - low).maxCoeff(&axis); + return axis; +} + +/*! \brief Throw on a segment center inside an unjoined segment. + \param views One view per final segment, reordered in place by the sweep. + \exception nec_exception* If two unconnected segments overlap. + + A center inside another segment's volume lies within that segment's own reach + of its center, so ordering the views along one axis bounds the window of + candidate pairs by the widest reach in the structure. +*/ +static void reject_overlapping_pairs(std::vector& views) +{ + nec_float reach_max = 0.0; + for (size_t i = 0; i < views.size(); i++) + reach_max = std::max(reach_max, + 0.5 * views[i].body.length() + views[i].body.get_radius()); + + Eigen::Index axis = widest_extent_axis(views); + std::stable_sort(views.begin(), views.end(), + [axis](const segment_view& lhs, const segment_view& rhs) { + return lhs.midpoint(axis) < rhs.midpoint(axis); + }); + + for (size_t anchor = 0; anchor < views.size(); anchor++) + { + for (size_t probe = anchor + 1; + probe < views.size() && + (views[probe].midpoint(axis) - views[anchor].midpoint(axis)) <= reach_max; + probe++) + { + // A recorded junction is a connection, not an overlap. + if (segments_joined(views[anchor], views[probe])) + continue; + + if (views[probe].body.intersect(views[anchor].midpoint)) + throw_segment_overlap(views[anchor], views[probe]); + + if (views[anchor].body.intersect(views[probe].midpoint)) + throw_segment_overlap(views[probe], views[anchor]); + } + } +} + +/* +There is no restriction on the angle between two wires, but accuracy will be lost if the center of a segment falls within the +volume of the wire the segment connects to. The risk of this reduces as the angle between wires approaches 180 degrees. + +Wires which intersect away from their ends are not connected, but errors will occur if one wire occupies the space of another one. For accuracy, separate wire centers by several radii of the largest wire. +*/ +void c_geometry::check_segment_intersections() +{ + if (false == _check_intersections) + return; + + if (n_segments < 2) + return; + + // NEC stores each final segment as a center, a length, and a direction, so + // the endpoints follow from a half-length step either side of the center. + std::vector views; + views.reserve(n_segments); + for (int64_t i = 0; i < n_segments; i++) + { + nec_3vector center(x[i], y[i], z[i]); + nec_3vector half_axis = nec_3vector(cab[i], sab[i], salp[i]) + * (0.5 * segment_length[i]); + + views.push_back(segment_view{ + i, + nec_wire(center - half_axis, center + half_axis, + segment_radius[i], segment_tags[i]), + center, + connected_segment_index(icon1[i]), + connected_segment_index(icon2[i])}); + } + + reject_overlapping_pairs(views); +} + void c_geometry::build_connections( int ignd ) { if ( ignd != 0) { @@ -1474,8 +1580,8 @@ void c_geometry::build_connections( int ignd ) nec_3vector v1(x[i], y[i], z[i]); nec_3vector v2(x2[i], y2[i], z2[i]); - nec_float slen = norm(v2 - v1) * SMIN; - + nec_float slen = nec_contact_threshold(v1, v2); + /* determine connection data for end 1 of segment. */ bool segment_on_ground = false; if ( ignd > 0) { @@ -1495,29 +1601,29 @@ void c_geometry::build_connections( int ignd ) if ( false == segment_on_ground ) { int ic= i; - nec_float sep=0.0; + bool contact_found = false; for (int64_t j = 1; j < n_segments; j++) { ic++; if ( ic >= n_segments) ic=0; - + nec_3vector vic(x[ic], y[ic], z[ic]); - sep = normL1(v1 - vic); - if ( sep <= slen) { + if ( nec_points_contact(v1, vic, slen) ) { icon1[i]= -(ic+1); + contact_found = true; break; } - + nec_3vector v2ic(x2[ic], y2[ic], z2[ic]); - sep = normL1(v1 - v2ic); - if ( sep <= slen) { + if ( nec_points_contact(v1, v2ic, slen) ) { icon1[i]= (ic+1); + contact_found = true; break; } - + } /* for( j = 1; j < n_segments; j++) */ - - if ( ((iz > 0) || (icon1[i] <= PCHCON)) && (sep > slen) ) + + if ( ((iz > 0) || (icon1[i] <= PCHCON)) && (false == contact_found) ) icon1[i]=0; } /* if ( ! jump ) */ @@ -1549,28 +1655,28 @@ void c_geometry::build_connections( int ignd ) v1 = nec_3vector(x[i], y[i], z[i]); v2 = nec_3vector(x2[i], y2[i], z2[i]); int ic= i; - nec_float sep=0.0; + bool contact_found = false; for (int64_t j = 1; j < n_segments; j++ ) { ic++; if ( ic >= n_segments) ic=0; - + nec_3vector vic(x[ic], y[ic], z[ic]); - sep = normL1(v2 - vic); - if (sep <= slen) { + if (nec_points_contact(v2, vic, slen)) { icon2[i]= (ic+1); + contact_found = true; break; } - + nec_3vector v2ic(x2[ic], y2[ic], z2[ic]); - sep = normL1(v2 - v2ic); - if (sep <= slen) { + if (nec_points_contact(v2, v2ic, slen)) { icon2[i]= -(ic+1); + contact_found = true; break; } } /* for( j = 1; j < n_segments; j++ ) */ - - if ( ((iz > 0) || (icon2[i] <= PCHCON)) && (sep > slen) ) + + if ( ((iz > 0) || (icon2[i] <= PCHCON)) && (false == contact_found) ) icon2[i]=0; } /* for( i = 0; i < n_segments; i++ ) */ diff --git a/src/c_geometry.h b/src/c_geometry.h index 545de87e..b2629b7f 100644 --- a/src/c_geometry.h +++ b/src/c_geometry.h @@ -27,7 +27,6 @@ class nec_context; #include "nec_output.h" -#include "nec_wire.h" @@ -208,6 +207,13 @@ class c_geometry void divide_patch( int nx ); void connect_segments( int ignd ); + + /*! \brief Reject geometry in which a segment center lies inside the volume of + * a segment it is not connected to. + * \exception nec_exception* If two unconnected segments overlap. + */ + void check_segment_intersections(); + void build_connections( int ignd ); void resolve_junctions(); @@ -232,8 +238,6 @@ class c_geometry nec_context* m_context; nec_output_file* m_output; bool _check_intersections; - - std::vector m_wires; void reflect_plane(int sym_plane, int& tag_increment); diff --git a/src/nec_context_tb.cpp b/src/nec_context_tb.cpp index 4966a915..3a2f3c34 100644 --- a/src/nec_context_tb.cpp +++ b/src/nec_context_tb.cpp @@ -211,6 +211,71 @@ 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( "Helix rejects invalid segment_count", "[helix]") { // Regression test for #48: helix with segment_count < 1 must throw // rather than silently returning with no wires. diff --git a/src/nec_wire.h b/src/nec_wire.h index 26523a38..1356ba6b 100644 --- a/src/nec_wire.h +++ b/src/nec_wire.h @@ -94,6 +94,16 @@ std::vector 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) From 05eb355b1affe448005f76b0ed4c8472ea98fc3a Mon Sep 17 00:00:00 2001 From: Eric Wheeler Date: Wed, 5 Aug 2026 23:37:13 -0700 Subject: [PATCH 2/3] geometry: class junction nodes geometrically and add an overlap warning channel The segment overlap check rejected legal NEC-2 junctions. It excluded a pair of segments from the volume test by reading icon1/icon2, which record at most one contact per segment end, so a node where six ends meet is recorded as a cycle rather than a clique and pairs not adjacent in that cycle looked unconnected. A ground-plane vertical on a mast reported GEOMETRY DATA ERROR -- SEGMENT #14 (TAG ID #1) MIDPOINT LIES WITHIN SEGMENT #66 (TAG ID #3) at a measured 0.01308 against a radius of 0.02500. Where an end meets the ground the record names the ground rather than a neighbour, so two legs meeting at z = 0 had no link at all. Classing segment ends into nodes from the final coordinates recovers the clique whatever order build_connections() linked them in, and the exclusion now applies to pairs sharing exactly one node: two segments holding both nodes are laid one along the other, not joined. - add nec_node_partition, a union-find over the 2*n_segments segment ends united on the larger of the two owning thresholds, the symmetric closure of the contact relation - replace segments_joined() and connected_segment_index() with shares_one_node(), so the check no longer reads the solver connection record - generalize widest_extent_axis() and extract sort_along_axis() over both the endpoint sweep and the segment sweep, which order positions by the same rule - record nec_overlap_finding per violation, carrying the distance and radius measured, so the diagnostic reports how far inside the conductor the center falls - add set_intersection_fatal(bool) and overlap_findings() beside the preserved set_intersection_check(bool), letting a caller read the measurement either policy produces - move the ROM2 step-size explanation latch from a file static to a context member, so a second context in one process still receives it - add four [segment_junction] regressions: the six-way node, the junction on the ground plane, the intrusion past a shared node, and the warning policy reading its findings back Signed-off-by: Eric Wheeler --- src/c_geometry.cpp | 329 ++++++++++++++++++++++++++++++----------- src/c_geometry.h | 38 ++++- src/nec_context.cpp | 3 +- src/nec_context.h | 5 + src/nec_context_tb.cpp | 67 +++++++++ 5 files changed, 355 insertions(+), 87 deletions(-) diff --git a/src/c_geometry.cpp b/src/c_geometry.cpp index 6e6d216c..51dcc290 100644 --- a/src/c_geometry.cpp +++ b/src/c_geometry.cpp @@ -26,19 +26,6 @@ #include #include -/*!\brief One final segment as the intersection sweep reads it. - - Carries the segment's own index so the sorted sweep still names it, its final - geometry as a wire, its center, and the segments its two ends are joined to. -*/ -struct segment_view { - int64_t index; - nec_wire body; - nec_3vector midpoint; - int64_t linked_start; - int64_t linked_end; -}; - /** * geometry_field_separator - Determine whether a character separates fields * @character: Character read from a geometry card @@ -69,6 +56,7 @@ c_geometry::c_geometry() maxcon = 0; _check_intersections = true; + _overlap_is_fatal = true; m_context = NULL; m_output = NULL; @@ -1390,71 +1378,84 @@ static bool nec_points_contact(const nec_3vector& a, const nec_3vector& b, } -/*! \brief Resolve one connection record to a segment index. - \param connection An icon1 or icon2 entry. - \return The zero-based index of the joined segment, or -1 when no segment is joined at that end. +/*! \brief One end of one segment, carried through the node partition sweep. - NEC encodes a free end as zero, a segment or ground connection as a signed - one-based segment number, and a patch connection above PCHCON. + The index numbers the ends of the whole structure, so end 1 of segment i is + 2*i and end 2 is 2*i+1. That numbering survives the sort, which reorders the + ends without renaming them. */ -static int64_t connected_segment_index(int32_t connection) -{ - int64_t magnitude = (connection < 0) - ? -static_cast(connection) - : static_cast(connection); - bool joins_segment = (0 != connection) && (magnitude <= PCHCON); - return joins_segment ? (magnitude - 1) : -1; -} +struct segment_endpoint { + nec_3vector position; + int64_t index; + nec_float threshold; +}; -/*! \brief Report whether NEC recorded these two segments as sharing an end. - \param a One segment view. - \param b The other segment view. - \return true if either segment records the other at either of its ends. +/*! \brief The two node classes touched by one segment's ends. +*/ +struct segment_nodes { + int64_t first; + int64_t second; +}; - build_connections() stops at the first contact found for each end, so a - junction of three or more segments is recorded as a chain rather than a - clique. Reading both ends of both segments covers every pair that chain - records. +/*! \brief The partition of every segment end into the node it belongs to. + + build_connections() records at most one contact per segment end, so a node + where three or more ends meet is recorded as a cycle rather than as a clique + and pairs that are not adjacent in that cycle look unconnected. Classing the + ends geometrically recovers the clique: every end in contact with any other + end of its class belongs to that class, whatever order the connection record + happened to link them in. */ -static bool segments_joined(const segment_view& a, const segment_view& b) +class nec_node_partition { - return (a.linked_start == b.index) || (a.linked_end == b.index) || - (b.linked_start == a.index) || (b.linked_end == a.index); -} +public: + nec_node_partition(std::vector& endpoints, int64_t segment_count); + + /*! \brief Report whether two segments meet at exactly one node. + \param segment_a One segment index. + \param segment_b The other segment index. + \return true when the two segments have one node class in common. + + One node in common is a junction: the currents of both segments are + continuous through it, which bounds the error a center inside the joined + wire introduces. Two segments holding both nodes in common are not joined + end to end but laid one along the other, so their volumes coincide over + their whole length and the pair carries no such continuity. + */ + bool shares_one_node(int64_t segment_a, int64_t segment_b) const; + +private: + int64_t root(int64_t endpoint_index) const; + void unite(int64_t lhs, int64_t rhs); -/*! \brief Report one segment center inside another segment. - \param inside The segment whose center falls within the other's volume. - \param container The segment whose volume contains that center. - \exception nec_exception* Always thrown, naming both segments. + mutable std::vector m_parent; + std::vector m_segment_nodes; +}; + +/*! \brief One final segment, carried through the overlap sweep. */ -static void throw_segment_overlap(const segment_view& inside, - const segment_view& container) -{ - nec_exception nex("GEOMETRY DATA ERROR -- SEGMENT #"); - nex.append(inside.index + 1); - nex.append(" (TAG ID #"); nex.append(inside.body.tag_id()); - nex.append(") MIDPOINT LIES WITHIN SEGMENT #"); - nex.append(container.index + 1); - nex.append(" (TAG ID #"); nex.append(container.body.tag_id()); - nex.append(")"); - throw nex; -} +struct segment_view { + int64_t index; + nec_wire body; + nec_3vector position; +}; -/*! \brief Find the axis over which the structure spreads furthest. - \param views One view per final segment. +/*! \brief Find the axis over which a set of positions spreads furthest. + \param items One view per segment, or one per segment end. \return The coordinate index to sweep along. Sweeping along this axis leaves the ordered window admitting the fewest candidate pairs. */ -static Eigen::Index widest_extent_axis(const std::vector& views) +template +static Eigen::Index widest_extent_axis(const std::vector& items) { - nec_3vector low = views[0].midpoint; + nec_3vector low = items[0].position; nec_3vector high = low; - for (size_t i = 1; i < views.size(); i++) + for (size_t i = 1; i < items.size(); i++) { - low = low.cwiseMin(views[i].midpoint); - high = high.cwiseMax(views[i].midpoint); + low = low.cwiseMin(items[i].position); + high = high.cwiseMax(items[i].position); } Eigen::Index axis = 0; @@ -1462,15 +1463,151 @@ static Eigen::Index widest_extent_axis(const std::vector& views) return axis; } -/*! \brief Throw on a segment center inside an unjoined segment. +/*! \brief Order a set of positions along one coordinate axis. + \param items The views to reorder in place. + \param axis The coordinate index chosen by widest_extent_axis(). +*/ +template +static void sort_along_axis(std::vector& items, Eigen::Index axis) +{ + std::stable_sort(items.begin(), items.end(), + [axis](const T& lhs, const T& rhs) { + return lhs.position(axis) < rhs.position(axis); + }); +} + +/*! \brief Class the ends of every final segment into the nodes they meet at. + \param endpoints Both ends of every segment, reordered by the sweep. + \param segment_count The number of final segments the ends belong to. +*/ +nec_node_partition::nec_node_partition(std::vector& endpoints, + int64_t segment_count) + : m_parent(endpoints.size()), m_segment_nodes(segment_count) +{ + for (size_t i = 0; i < m_parent.size(); i++) + m_parent[i] = int64_t(i); + + nec_float threshold_max = 0.0; + for (size_t i = 0; i < endpoints.size(); i++) + threshold_max = std::max(threshold_max, endpoints[i].threshold); + + Eigen::Index axis = widest_extent_axis(endpoints); + sort_along_axis(endpoints, axis); + + // The connection rule measures a gap against the threshold of the segment + // whose end is under test, which makes it asymmetric. A node is a property + // shared by both ends, so the partition unites on the larger of the two + // thresholds: the symmetric closure of that relation. + // + // An axis gap already exceeds the L1 separation, so once the gap passes the + // widest threshold in the structure no later end in the order can touch this + // one. + for (size_t anchor = 0; anchor < endpoints.size(); anchor++) + { + for (size_t probe = anchor + 1; + probe < endpoints.size() && + (endpoints[probe].position(axis) - endpoints[anchor].position(axis)) <= threshold_max; + probe++) + { + nec_float threshold = std::max(endpoints[anchor].threshold, + endpoints[probe].threshold); + if (nec_points_contact(endpoints[anchor].position, endpoints[probe].position, + threshold)) + unite(endpoints[anchor].index, endpoints[probe].index); + } + } + + for (int64_t i = 0; i < segment_count; i++) + m_segment_nodes[i] = segment_nodes{root(2*i), root(2*i + 1)}; +} + +/*! \brief The class representative of one segment end. + \param endpoint_index The structure-wide end number. + \return The end that names the class, compressing the path walked to reach it. +*/ +int64_t nec_node_partition::root(int64_t endpoint_index) const +{ + int64_t node = endpoint_index; + while (m_parent[node] != node) + { + m_parent[node] = m_parent[m_parent[node]]; + node = m_parent[node]; + } + return node; +} + +/*! \brief Merge the classes of two segment ends. + \param lhs One end. + \param rhs The other end. + + Two ends already in one class name the same representative, and naming it + again leaves the class unchanged. +*/ +void nec_node_partition::unite(int64_t lhs, int64_t rhs) +{ + m_parent[root(rhs)] = root(lhs); +} + +bool nec_node_partition::shares_one_node(int64_t segment_a, int64_t segment_b) const +{ + const segment_nodes& a = m_segment_nodes[segment_a]; + const segment_nodes& b = m_segment_nodes[segment_b]; + + int shared = int((a.first == b.first) || (a.first == b.second)) + + int((a.second == b.first) || (a.second == b.second)); + + return 1 == shared; +} + +/*! \brief Render one overlap finding as a diagnostic line. + \param prefix The severity the caller applies to the finding. + \param finding The measurement to render. + \return The rendered line, without a trailing newline. +*/ +static std::string nec_overlap_message(const char* prefix, + const nec_overlap_finding& finding) +{ + char buffer[256]; + snprintf(buffer, sizeof(buffer), + "%s -- SEGMENT #%lld (TAG ID #%d) MIDPOINT LIES WITHIN SEGMENT #%lld" + " (TAG ID #%d) AT DISTANCE %.5f RADIUS %.5f", + prefix, + (long long)(finding.inside_index + 1), finding.inside_tag, + (long long)(finding.container_index + 1), finding.container_tag, + finding.distance, finding.container_radius); + return std::string(buffer); +} + +/*! \brief Record one segment center falling inside another segment's volume. + \param inside The segment whose center is measured. + \param container The segment whose volume the center is measured against. + \param findings Appended with the measurement when the center lies inside. +*/ +static void append_overlap(const segment_view& inside, const segment_view& container, + std::vector& findings) +{ + nec_float distance = container.body.axis_distance(inside.position); + if (distance > container.body.get_radius()) + return; + + findings.push_back(nec_overlap_finding{ + inside.index, container.index, + inside.body.tag_id(), container.body.tag_id(), + distance, container.body.get_radius()}); +} + +/*! \brief Collect every segment center that lies inside an unjoined segment. \param views One view per final segment, reordered in place by the sweep. - \exception nec_exception* If two unconnected segments overlap. + \param nodes The node classing of the segment ends. + \param findings Appended with one record per violation. A center inside another segment's volume lies within that segment's own reach of its center, so ordering the views along one axis bounds the window of candidate pairs by the widest reach in the structure. */ -static void reject_overlapping_pairs(std::vector& views) +static void collect_overlapping_pairs(std::vector& views, + const nec_node_partition& nodes, + std::vector& findings) { nec_float reach_max = 0.0; for (size_t i = 0; i < views.size(); i++) @@ -1478,27 +1615,21 @@ static void reject_overlapping_pairs(std::vector& views) 0.5 * views[i].body.length() + views[i].body.get_radius()); Eigen::Index axis = widest_extent_axis(views); - std::stable_sort(views.begin(), views.end(), - [axis](const segment_view& lhs, const segment_view& rhs) { - return lhs.midpoint(axis) < rhs.midpoint(axis); - }); + sort_along_axis(views, axis); for (size_t anchor = 0; anchor < views.size(); anchor++) { for (size_t probe = anchor + 1; probe < views.size() && - (views[probe].midpoint(axis) - views[anchor].midpoint(axis)) <= reach_max; + (views[probe].position(axis) - views[anchor].position(axis)) <= reach_max; probe++) { - // A recorded junction is a connection, not an overlap. - if (segments_joined(views[anchor], views[probe])) + // Segments meeting at one node are joined, not overlapping. + if (nodes.shares_one_node(views[anchor].index, views[probe].index)) continue; - if (views[probe].body.intersect(views[anchor].midpoint)) - throw_segment_overlap(views[anchor], views[probe]); - - if (views[anchor].body.intersect(views[probe].midpoint)) - throw_segment_overlap(views[probe], views[anchor]); + append_overlap(views[anchor], views[probe], findings); + append_overlap(views[probe], views[anchor], findings); } } } @@ -1511,32 +1642,64 @@ Wires which intersect away from their ends are not connected, but errors will oc */ void c_geometry::check_segment_intersections() { + m_overlap_findings.clear(); + if (false == _check_intersections) return; if (n_segments < 2) return; + collect_overlap_findings(); + + // A fatal policy rejects the geometry on its first finding; a warning policy + // reports every finding and leaves the geometry in place. + if (_overlap_is_fatal && false == m_overlap_findings.empty()) + throw nec_exception(nec_overlap_message("GEOMETRY DATA ERROR", + m_overlap_findings.front()).c_str()); + else + { + // A fatal policy reaches this loop only with nothing to report. + for (size_t i = 0; i < m_overlap_findings.size(); i++) + m_output->nec_printf("\n%s", + nec_overlap_message("GEOMETRY WARNING", m_overlap_findings[i]).c_str()); + } +} + +/*! \brief Measure every segment center against the volume of every segment it + does not meet at a node. + + Appends one finding per violation to m_overlap_findings, which the caller has + already emptied. +*/ +void c_geometry::collect_overlap_findings() +{ // NEC stores each final segment as a center, a length, and a direction, so // the endpoints follow from a half-length step either side of the center. std::vector views; - views.reserve(n_segments); + views.reserve(size_t(n_segments)); + std::vector endpoints; + endpoints.reserve(size_t(2 * n_segments)); for (int64_t i = 0; i < n_segments; i++) { nec_3vector center(x[i], y[i], z[i]); nec_3vector half_axis = nec_3vector(cab[i], sab[i], salp[i]) * (0.5 * segment_length[i]); + nec_3vector end_one = center - half_axis; + nec_3vector end_two = center + half_axis; views.push_back(segment_view{ i, - nec_wire(center - half_axis, center + half_axis, - segment_radius[i], segment_tags[i]), - center, - connected_segment_index(icon1[i]), - connected_segment_index(icon2[i])}); + nec_wire(end_one, end_two, segment_radius[i], segment_tags[i]), + center}); + + nec_float threshold = nec_contact_threshold(end_one, end_two); + endpoints.push_back(segment_endpoint{end_one, 2*i, threshold}); + endpoints.push_back(segment_endpoint{end_two, 2*i + 1, threshold}); } - reject_overlapping_pairs(views); + nec_node_partition nodes(endpoints, n_segments); + collect_overlapping_pairs(views, nodes, m_overlap_findings); } void c_geometry::build_connections( int ignd ) diff --git a/src/c_geometry.h b/src/c_geometry.h index b2629b7f..b7ce7d43 100644 --- a/src/c_geometry.h +++ b/src/c_geometry.h @@ -29,6 +29,20 @@ class nec_context; #include "nec_output.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 @@ -133,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& overlap_findings() const { return m_overlap_findings; } + /*! \brief Geometry is complete * \exception nec_exception* If there is an error with the geometry. */ @@ -208,12 +234,16 @@ class c_geometry void connect_segments( int ignd ); - /*! \brief Reject geometry in which a segment center lies inside the volume of - * a segment it is not connected to. - * \exception nec_exception* If two unconnected segments overlap. + /*! \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(); @@ -238,6 +268,8 @@ class c_geometry nec_context* m_context; nec_output_file* m_output; bool _check_intersections; + bool _overlap_is_fatal; + std::vector m_overlap_findings; void reflect_plane(int sym_plane, int& tag_increment); diff --git a/src/nec_context.cpp b/src/nec_context.cpp index df5faf64..52ef7ae7 100644 --- a/src/nec_context.cpp +++ b/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; diff --git a/src/nec_context.h b/src/nec_context.h index 6aeec19c..a4a8c757 100644 --- a/src/nec_context.h +++ b/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); diff --git a/src/nec_context_tb.cpp b/src/nec_context_tb.cpp index 3a2f3c34..3cbe3a1c 100644 --- a/src/nec_context_tb.cpp +++ b/src/nec_context_tb.cpp @@ -276,6 +276,73 @@ TEST_CASE( "Wire penetrating another away from any node is rejected", "[segment_ 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. From 6910f284e98d71852c4205f27c8a056d986b6488 Mon Sep 17 00:00:00 2001 From: Eric Wheeler Date: Wed, 5 Aug 2026 23:37:38 -0700 Subject: [PATCH 3/3] geometry: add a -w option selecting the overlap warning policy The overlap policy introduced with the node classing was reachable only from the C++ library, so a command-line user with a model that trips the volume test had no way to see the measurement and continue. This exposes that policy as a command-line option, leaving the default unchanged: without -w the first finding still rejects the geometry. - accept -w in the option string and clear the geometry fatal flag when it is given, so each finding is reported after the segment table and the solve proceeds - name the option in the usage text, so its effect on a rejected model is discoverable Signed-off-by: Eric Wheeler --- src/misc.cpp | 1 + src/nec2cpp.cpp | 6 +++++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/misc.cpp b/src/misc.cpp index f8208edd..b56ad0c5 100644 --- a/src/misc.cpp +++ b/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"; diff --git a/src/nec2cpp.cpp b/src/nec2cpp.cpp index f1a4cf9b..707121a8 100644 --- a/src/nec2cpp.cpp +++ b/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);