Skip to content

[WIP] [WIP] Address feedback to apply coefficient rotation change in computeAccel() calls #221

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
michael-petersen merged 5 commits into from
May 29, 2026
Merged

[WIP] [WIP] Address feedback to apply coefficient rotation change in computeAccel() calls #221

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9256441
Initial plan
Copilot May 27, 2026
b4470e1
Apply coefrot in computeAccel implementations
Copilot May 27, 2026
19f61d3
Fix FlatDisk off-grid component order in rotated accel
Copilot May 27, 2026
407afbf
Merge branch 'devel' into copilot/wip-address-feedback-apply-coeffici…
michael-petersen May 28, 2026
a935940
Fix Cube::computeAccel redeclaration from coefficient rotation change
Copilot May 28, 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
94 changes: 47 additions & 47 deletions expui/BiorthBasis.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -818,19 +818,10 @@ namespace BasisClasses
void Spherical::computeAccel(double x, double y, double z,
Eigen::Ref<Eigen::Vector3d> acc)
{
struct AccelFrameGuard {
Eigen::Ref<Eigen::Vector3d> acc;
const decltype(coefrot)& rot;

~AccelFrameGuard()
{
acc = rot.transpose() * acc;
}
} rotate_acc_to_caller{acc, coefrot};

// Shift to center and rotate into the coefficient frame
Eigen::Vector3d pos(x - coefctr(0), y - coefctr(1), z - coefctr(2));
pos = coefrot * pos;
// Shift and rotate to expansion frame
Eigen::Vector3d pos {x, y, z};
pos = coefrot * (pos - coefctr);

x = pos(0);
y = pos(1);
z = pos(2);
Expand Down Expand Up @@ -937,9 +928,9 @@ namespace BasisClasses
double tpoty = (potr - pott*costh/r)*y/r + potp*x/R2;
double tpotz = potr*costh + pott*sinth*sinth/r;

// Return force not potential gradient
// Return force not potential gradient and rotate to caller frame
//
acc << tpotx, tpoty, tpotz;
acc = coefrot.transpose() * Eigen::Vector3d(tpotx, tpoty, tpotz);
}


Expand Down Expand Up @@ -1825,10 +1816,12 @@ namespace BasisClasses
void Cylindrical::computeAccel(double x, double y, double z,
Eigen::Ref<Eigen::Vector3d> acc)
{
// Shift to center
x -= coefctr(0);
y -= coefctr(1);
z -= coefctr(2);
// Shift and rotate to expansion frame
Eigen::Vector3d pos {x, y, z};
pos = coefrot * (pos - coefctr);
x = pos(0);
y = pos(1);
z = pos(2);

double R = sqrt(x*x + y*y);
double phi = atan2(y, x);
Expand All @@ -1846,8 +1839,8 @@ namespace BasisClasses
tpoty = tpotR*y/R + tpotp*x/R;
}

// Apply G to forces on return
acc << tpotx*G, tpoty*G, tpotz*G;
// Apply G and rotate forces on return
acc = coefrot.transpose() * Eigen::Vector3d(tpotx*G, tpoty*G, tpotz*G);
}

// Evaluate in cylindrical coordinates
Expand Down Expand Up @@ -2511,10 +2504,12 @@ namespace BasisClasses
void FlatDisk::computeAccel(double x, double y, double z,
Eigen::Ref<Eigen::Vector3d> acc)
{
// Shift to center
x -= coefctr(0);
y -= coefctr(1);
z -= coefctr(2);
// Shift and rotate to expansion frame
Eigen::Vector3d pos {x, y, z};
pos = coefrot * (pos - coefctr);
x = pos(0);
y = pos(1);
z = pos(2);

// Get thread id
int tid = omp_get_thread_num();
Expand All @@ -2539,7 +2534,7 @@ namespace BasisClasses
rpot = -G*totalMass*R/(r*r2 + 10.0*std::numeric_limits<double>::min());
zpot = -G*totalMass*z/(r*r2 + 10.0*std::numeric_limits<double>::min());

acc << rpot, zpot, ppot;
acc = coefrot.transpose() * Eigen::Vector3d(rpot, ppot, zpot);
}

// Get the basis fields
Expand Down Expand Up @@ -2607,7 +2602,7 @@ namespace BasisClasses
poty = rpot*y/R + ppot*x/R;
}

acc << potx, poty, zpot;
acc = coefrot.transpose() * Eigen::Vector3d(potx, poty, zpot);
}


Expand Down Expand Up @@ -3309,10 +3304,12 @@ namespace BasisClasses
void CBDisk::computeAccel(double x, double y, double z,
Eigen::Ref<Eigen::Vector3d> acc)
{
// Shift to center
x -= coefctr(0);
y -= coefctr(1);
z -= coefctr(2);
// Shift and rotate to expansion frame
Eigen::Vector3d pos {x, y, z};
pos = coefrot * (pos - coefctr);
x = pos(0);
y = pos(1);
z = pos(2);

// Get thread id
int tid = omp_get_thread_num();
Expand Down Expand Up @@ -3381,7 +3378,7 @@ namespace BasisClasses
poty = rpot*y/R + ppot*x/R;
}

acc << potx, poty, zpot;
acc = coefrot.transpose() * Eigen::Vector3d(potx, poty, zpot);
}


Expand Down Expand Up @@ -3804,10 +3801,12 @@ namespace BasisClasses
void Slab::computeAccel(double x, double y, double z,
Eigen::Ref<Eigen::Vector3d> acc)
{
// Shift to center
x -= coefctr(0);
y -= coefctr(1);
z -= coefctr(2);
// Shift and rotate to expansion frame
Eigen::Vector3d pos {x, y, z};
pos = coefrot * (pos - coefctr);
x = pos(0);
y = pos(1);
z = pos(2);

// Loop indices
//
Expand Down Expand Up @@ -3880,8 +3879,9 @@ namespace BasisClasses
}
}

// Apply G to forces on return
acc << G*accx.real(), G*accy.real(), G*accz.real();
// Apply G to forces on return and rotate to caller frame
acc = coefrot.transpose() *
Eigen::Vector3d(G*accx.real(), G*accy.real(), G*accz.real());
}


Expand Down Expand Up @@ -4381,22 +4381,22 @@ namespace BasisClasses
void Cube::computeAccel(double x, double y, double z,
Eigen::Ref<Eigen::Vector3d> acc)
{
// Shift to center
x -= coefctr(0);
y -= coefctr(1);
z -= coefctr(2);
// Shift and rotate to expansion frame
Eigen::Vector3d pos {x, y, z};
pos = coefrot * (pos - coefctr);
x = pos(0);
y = pos(1);
z = pos(2);

// Get thread id
int tid = omp_get_thread_num();

// Position vector
Eigen::Vector3d pos {x, y, z};

// Get the basis fields
auto frc = ortho->get_force(expcoef, pos);

// Apply G to forces on return
acc << -G*frc(0).real(), -G*frc(1).real(), -G*frc(2).real();
// Apply G to forces on return and rotate to caller frame
acc = coefrot.transpose() *
Eigen::Vector3d(-G*frc(0).real(), -G*frc(1).real(), -G*frc(2).real());
}

std::vector<double> Cube::cyl_eval(double R, double z, double phi)
Expand Down
Loading