Skip to content

symmetry: fix gx solve and three-plane reflection #124

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
tmolteno merged 3 commits into from
Aug 2, 2026
+44 −62
Merged

symmetry: fix gx solve and three-plane reflection #124

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
af8b573
fix: seed scm[0] before symmetry-mode sum overwrites b[i]
Aug 1, 2026
33f93d8
fix: correct plane-symmetry matrix fill guard in fblock()
Aug 1, 2026
7475b20
fix: reflect three symmetry planes in sequential Z, Y, X order
Aug 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
75 changes: 23 additions & 52 deletions src/c_geometry.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1152,29 +1152,28 @@ void c_geometry::reflect_plane( int sym_plane, int& tag_increment )

/* --- planar symmetry --- */
/* sym_plane bits: 0 = negate X, 1 = negate Y, 2 = negate Z */
int mask = sym_plane;
bool neg_x = mask & 1;
bool neg_y = mask & 2;
bool neg_z = mask & 4;
int num_axes = (neg_x?1:0) + (neg_y?1:0) + (neg_z?1:0);
int num_copies = 1 << num_axes; /* 2 for single axis, 4 for two axes */

int64_t orig_n = n_segments;
int64_t orig_m = m;
int itx = tag_increment;

/* Axis priority: Z (bit 2) > Y (bit 1) > X (bit 0).
For two-axis combos, copy order is: original, axis1-neg, axis2-neg, both-neg.
Tag increments: +0, +itx, +2*itx, +3*itx. */
int axis1_mask = 0, axis2_mask = 0;
if ( neg_z ) { axis1_mask = 4; axis2_mask = neg_y ? 2 : (neg_x ? 1 : 0); }
else if ( neg_y ) { axis1_mask = 2; axis2_mask = neg_x ? 1 : 0; }
else if ( neg_x ) { axis1_mask = 1; }
/* Build copy masks in the Z, Y, X order used by sequential NEC-2
reflection passes. Each pass appends a reflected copy of every existing
combination. */
int copy_mask[8] = { 0 };
int num_copies = 1;

for ( int axis_mask = 4; axis_mask > 0; axis_mask >>= 1 )
{
if ( 0 == (sym_plane & axis_mask) )
continue;

int tag_inc[4] = { 0, itx, 2*itx, 3*itx };
for ( int copy = 0; copy < num_copies; copy++ )
copy_mask[num_copies + copy] = copy_mask[copy] | axis_mask;

num_copies *= 2;
}

/* Final tag_increment: X never doubles, Z and Y always double */
tag_increment = itx << ( (neg_z ? 1 : 0) + (neg_y ? 1 : 0) );
tag_increment = itx * num_copies;

/* --- SEGMENTS --- */
if ( orig_n > 0 )
Expand All @@ -1196,23 +1195,9 @@ void c_geometry::reflect_plane( int sym_plane, int& tag_increment )
for( int copy = 1; copy < num_copies; copy++ )
{
/* Determine which axes to negate for this copy */
bool flip_x = false, flip_y = false, flip_z = false;
if ( copy == 1 )
{
if ( axis1_mask == 4 ) flip_z = true;
else if ( axis1_mask == 2 ) flip_y = true;
else if ( axis1_mask == 1 ) flip_x = true;
}
else if ( copy == 2 )
{
if ( axis2_mask == 4 ) flip_z = true;
else if ( axis2_mask == 2 ) flip_y = true;
else if ( axis2_mask == 1 ) flip_x = true;
}
else if ( copy == 3 )
{
flip_x = neg_x; flip_y = neg_y; flip_z = neg_z;
}
bool flip_x = (0 != (copy_mask[copy] & 1));
bool flip_y = (0 != (copy_mask[copy] & 2));
bool flip_z = (0 != (copy_mask[copy] & 4));

int64_t base = copy * orig_n;

Expand Down Expand Up @@ -1269,7 +1254,7 @@ void c_geometry::reflect_plane( int sym_plane, int& tag_increment )
if ( itagi == 0 )
segment_tags[nx] = 0;
if ( itagi != 0 )
segment_tags[nx] = itagi + tag_inc[copy];
segment_tags[nx] = itagi + copy * itx;

segment_radius[nx] = segment_radius[i];
}
Expand Down Expand Up @@ -1298,23 +1283,9 @@ void c_geometry::reflect_plane( int sym_plane, int& tag_increment )

for( int copy = 1; copy < num_copies; copy++ )
{
bool flip_x = false, flip_y = false, flip_z = false;
if ( copy == 1 )
{
if ( axis1_mask == 4 ) flip_z = true;
else if ( axis1_mask == 2 ) flip_y = true;
else if ( axis1_mask == 1 ) flip_x = true;
}
else if ( copy == 2 )
{
if ( axis2_mask == 4 ) flip_z = true;
else if ( axis2_mask == 2 ) flip_y = true;
else if ( axis2_mask == 1 ) flip_x = true;
}
else if ( copy == 3 )
{
flip_x = neg_x; flip_y = neg_y; flip_z = neg_z;
}
bool flip_x = (0 != (copy_mask[copy] & 1));
bool flip_y = (0 != (copy_mask[copy] & 2));
bool flip_z = (0 != (copy_mask[copy] & 4));

int neg_count = (flip_x?1:0) + (flip_y?1:0) + (flip_z?1:0);
bool neg_psalp = (neg_count % 2 == 1);
Expand Down
12 changes: 10 additions & 2 deletions src/matrix_algebra.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -398,7 +398,11 @@ void solves(complex_array& a, int_array& ip, complex_array& b, int64_t neq,

/* transform matrix eq. rhs vector according to symmetry modes */
for (int64_t i = 0; i < npeq; i++ ) {
nec_complex sum_normal(b[i+column_offset]);
/* Each image row below is seeded from this row's untransformed
value, so hold it in scm[0] before the normalized sum
overwrites b[i]. */
scm[0] = b[i+column_offset];
nec_complex sum_normal(scm[0]);
for (int64_t k = 1; k < nop; k++ ) {
int64_t ia= i+ k* npeq;
scm[k]= b[ia+column_offset];
Expand Down Expand Up @@ -441,7 +445,11 @@ void solves(complex_array& a, int_array& ip, complex_array& b, int64_t neq,
for (int64_t ic = 0; ic < nrh; ic++ ) {
int64_t column_offset = ic*neq;
for (int64_t i = 0; i < npeq; i++ ) {
nec_complex sum_normal(b[i+column_offset]);
/* Each image row below is seeded from this row's untransformed
value, so hold it in scm[0] before the summed value overwrites
b[i]. */
scm[0] = b[i+column_offset];
nec_complex sum_normal(scm[0]);
for (int64_t k = 1; k < nop; k++ ) {
int64_t ia= i+ k* npeq;
scm[k]= b[ia+column_offset];
Expand Down
19 changes: 11 additions & 8 deletions src/nec_context.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6439,14 +6439,17 @@ void nec_context::fblock( int nrow, int ncol, int64_t imax, int ipsym ) {
int kk=1;
symmetry_array[0]=cplx_10();

if ((2 == nop) || (4 == nop) || (8 == nop))
return;
int ka = nop / 2;

// int k_power = 2;
// for( ka = 1; k_power != nop; ka++ )
// k_power *= 2;

/* Plane symmetry admits two, four, or eight sections; the doubling below
cannot span the matrix for any other count. */
if ((2 != nop) && (4 != nop) && (8 != nop))
nec_stop("SYMMETRY ERROR - NOP: %d", nop );

/* Each pass doubles the width of the filled block, so the pass count is
log2(nop). */
int ka = 1;
for (int k_power = 2; k_power != nop; k_power *= 2)
ka++;

for(int k = 0; k < ka; k++ ) {
for(int i = 0; i < kk; i++ ) {
for(int j = 0; j < kk; j++ ) {
Expand Down
Loading