Skip to content

Commit 2b31925

Browse files
michael-petersenmichael-petersen
authored
Merge pull request #153 from EXP-code/RotationMatrix
Implement rotation matrices for frame transformation
2 parents c9aff75 + 8c48318 commit 2b31925

18 files changed

Lines changed: 440 additions & 156 deletions

‎expui/BasisFactory.H‎

Lines changed: 19 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,8 @@ namespace BasisClasses
2424
{
2525
using RowMatrixXd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
2626

27+
using RowMatrix3d = Eigen::Matrix<double, 3, 3, Eigen::RowMajor>;
28+
2729
//! Callback function signature for selection particles to
2830
//! accumulate
2931
using Callback =
@@ -66,7 +68,10 @@ namespace BasisClasses
6668
unsigned long coefindx;
6769

6870
//! The expansion center
69-
std::vector<double> coefctr;
71+
Eigen::Vector3d coefctr = Eigen::Vector3d::Zero();
72+
73+
//! Rotation matrix
74+
RowMatrix3d coefrot = Eigen::Matrix3d::Identity();
7075

7176
//! Contains contructor and BFE parameter database
7277
YAML::Node node, conf;
@@ -151,8 +156,10 @@ namespace BasisClasses
151156
virtual ~Basis(void) {}
152157

153158
//! Set the expansion center
154-
void setCenter(std::vector<double> center)
155-
{ coefctr = center; }
159+
void setCenter(Eigen::Vector3d center) { coefctr = center; }
160+
161+
//! Set the rotation matrix
162+
void setRotation(RowMatrix3d rot) { coefrot = rot; }
156163

157164
//! Evaluate basis in desired coordinates
158165
virtual std::vector<double>
@@ -213,11 +220,12 @@ namespace BasisClasses
213220

214221
//! Generate coeffients from a particle reader
215222
virtual CoefClasses::CoefStrPtr
216-
createFromReader(PR::PRptr reader, std::vector<double> ctr) = 0;
223+
createFromReader(PR::PRptr reader, Eigen::Vector3d ctr,
224+
RowMatrix3d rot) = 0;
217225

218226
//! Generate coefficients from a phase-space table
219227
virtual void
220-
initFromArray(std::vector<double> ctr) = 0;
228+
initFromArray(Eigen::Vector3d ctr, RowMatrix3d rot) = 0;
221229

222230
//! Accumulate coefficient contributions from arrays
223231
virtual void
@@ -228,7 +236,8 @@ namespace BasisClasses
228236
//! for the expansion
229237
CoefClasses::CoefStrPtr createFromArray
230238
(Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0,
231-
std::vector<double> center={0.0, 0.0, 0.0},
239+
Eigen::Vector3d center=Eigen::Vector3d::Zero(),
240+
RowMatrix3d rot=RowMatrix3d::Identity(),
232241
bool roundrobin=true, bool posvelrows=false);
233242

234243
//! Create and the coefficients from the array accumulation with the
@@ -275,7 +284,10 @@ namespace BasisClasses
275284
{ return getFieldLabels(coordinates); }
276285

277286
//! Get the basis expansion center
278-
std::vector<double> getCenter() { return coefctr; }
287+
Eigen::Vector3d getCenter() { return coefctr; }
288+
289+
//! Get the basis expansion center
290+
RowMatrix3d getRotation() { return coefrot; }
279291
};
280292

281293
using BasisPtr = std::shared_ptr<Basis>;

‎expui/BasisFactory.cc‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -275,10 +275,10 @@ namespace BasisClasses
275275
// Generate coefficients from a phase-space table
276276
//
277277
CoefClasses::CoefStrPtr Basis::createFromArray
278-
(Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector<double> ctr,
279-
bool roundrobin, bool posvelrows)
278+
(Eigen::VectorXd& m, RowMatrixXd& p, double time, Eigen::Vector3d ctr,
279+
RowMatrix3d rot, bool roundrobin, bool posvelrows)
280280
{
281-
initFromArray(ctr);
281+
initFromArray(ctr, rot);
282282
addFromArray(m, p, roundrobin, posvelrows);
283283
return makeFromArray(time);
284284
}

‎expui/BiorthBasis.H‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -112,13 +112,15 @@ namespace BasisClasses
112112
//! Generate coeffients from a particle reader and optional center
113113
//! location for the expansion
114114
CoefClasses::CoefStrPtr createFromReader
115-
(PR::PRptr reader, std::vector<double> center={0.0, 0.0, 0.0});
115+
(PR::PRptr reader, Eigen::Vector3d center=Eigen::Vector3d::Zero(),
116+
RowMatrix3d rot=RowMatrix3d::Identity());
116117

117118
//! Generate coeffients from an array and optional center location
118119
//! for the expansion
119120
CoefClasses::CoefStrPtr createFromArray
120121
(Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0,
121-
std::vector<double> center={0.0, 0.0, 0.0},
122+
Eigen::Vector3d center=Eigen::Vector3d::Zero(),
123+
RowMatrix3d rot=RowMatrix3d::Identity(),
122124
bool roundrobin=true, bool posvelrows=false);
123125

124126
//! Generate coeffients from an array and optional center location
@@ -128,7 +130,8 @@ namespace BasisClasses
128130
//! Initialize accumulating coefficients from arrays with an optional
129131
//! center vector. This is called once to initialize the accumulation.
130132
void initFromArray
131-
(std::vector<double> center={0.0, 0.0, 0.0});
133+
(Eigen::Vector3d center=Eigen::Vector3d::Zero(),
134+
RowMatrix3d rot=RowMatrix3d::Identity());
132135

133136
//! Initialize accumulating coefficients from arrays. This is
134137
//! called once to initialize the accumulation.

‎expui/BiorthBasis.cc‎

Lines changed: 94 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -3956,7 +3956,7 @@ namespace BasisClasses
39563956

39573957
// Generate coeffients from a particle reader
39583958
CoefClasses::CoefStrPtr BiorthBasis::createFromReader
3959-
(PR::PRptr reader, std::vector<double> ctr)
3959+
(PR::PRptr reader, Eigen::Vector3d ctr, RowMatrix3d rot)
39603960
{
39613961
CoefClasses::CoefStrPtr coef;
39623962

@@ -3975,44 +3975,55 @@ namespace BasisClasses
39753975
throw std::runtime_error(sout.str());
39763976
}
39773977

3978-
// Is center non-zero?
3978+
// Add the expansion center metadata and register for this instance
39793979
//
3980-
bool addCenter = false;
3981-
for (auto v : ctr) {
3982-
if (v != 0.0) addCenter = true;
3983-
}
3980+
coefctr = ctr;
3981+
coef->ctr = ctr;
39843982

3985-
// Add the expansion center metadata
3983+
// Add the rotation matrix metadata and register for this instance
39863984
//
3987-
if (addCenter) coef->ctr = ctr;
3985+
coefrot = rot;
3986+
coef->rot = rot;
39883987

3989-
std::vector<double> pp(3), vv(3);
3988+
std::vector<double> p1(3), v1(3);
3989+
3990+
// Map the vector rather than copy
3991+
//
3992+
Eigen::Map<Eigen::Vector3d> pp(p1.data(), 3), vv(v1.data(), 3);
3993+
vv.setZero();
39903994

39913995
reset_coefs();
39923996
for (auto p=reader->firstParticle(); p!=0; p=reader->nextParticle()) {
39933997

39943998
bool use = false;
39953999

4000+
// Translate and rotate the position vector
4001+
//
4002+
for (int k=0; k<3; k++) p1[k] = p->pos[k];
4003+
pp = coefrot * (pp - coefctr);
4004+
39964005
if (ftor) {
3997-
pp.assign(p->pos, p->pos+3);
3998-
vv.assign(p->vel, p->vel+3);
3999-
use = ftor(p->mass, pp, vv, p->indx);
4006+
// Rotate the velocity vector
4007+
//
4008+
for (int k=0; k<3; k++) v1[k] = p->vel[k];
4009+
vv = coefrot * vv;
4010+
4011+
use = ftor(p->mass, p1, v1, p->indx);
40004012
} else {
40014013
use = true;
40024014
}
40034015

4004-
if (use) accumulate(p->pos[0]-ctr[0],
4005-
p->pos[1]-ctr[1],
4006-
p->pos[2]-ctr[2],
4007-
p->mass);
4016+
if (use) {
4017+
accumulate(pp(0), pp(1), pp(2), p->mass);
4018+
}
40084019
}
40094020
make_coefs();
40104021
load_coefs(coef, reader->CurrentTime());
40114022
return coef;
40124023
}
40134024

40144025
// Generate coefficients from a phase-space table
4015-
void BiorthBasis::initFromArray(std::vector<double> ctr)
4026+
void BiorthBasis::initFromArray(Eigen::Vector3d ctr, RowMatrix3d rot)
40164027
{
40174028
if (name.compare("sphereSL") == 0)
40184029
coefret = std::make_shared<CoefClasses::SphStruct>();
@@ -4027,20 +4038,14 @@ namespace BasisClasses
40274038
throw std::runtime_error(sout.str());
40284039
}
40294040

4030-
// Is center non-zero?
4031-
//
4032-
bool addCenter = false;
4033-
for (auto v : ctr) {
4034-
if (v != 0.0) addCenter = true;
4035-
}
4036-
4037-
// Add the expansion center metadata
4038-
//
4039-
if (addCenter) coefret->ctr = ctr;
4040-
4041-
// Register the center
4041+
// Add the expansion center metadata and register
40424042
//
40434043
coefctr = ctr;
4044+
coefret->ctr = ctr;
4045+
4046+
// Add the rotation metadata and register
4047+
coefrot = rot;
4048+
coefret->rot = rot;
40444049

40454050
// Clean up for accumulation
40464051
//
@@ -4069,6 +4074,7 @@ namespace BasisClasses
40694074
int cols = p.cols();
40704075

40714076
bool ambiguous = false;
4077+
bool haveVel = false;
40724078

40734079
if (cols==3 or cols==6) {
40744080
if (rows != 3 and rows != 6) PosVelRows = false;
@@ -4089,7 +4095,11 @@ namespace BasisClasses
40894095
<< "if this assumption is wrong." << std::endl;
40904096
}
40914097

4092-
std::vector<double> p1(3), v1(3, 0);
4098+
// Map the vector rather than copy
4099+
//
4100+
std::vector<double> p1(3), v1(3);
4101+
Eigen::Map<Eigen::Vector3d> pp(p1.data(), 3), vv(v1.data(), 3);
4102+
vv.setZero();
40934103

40944104
if (PosVelRows) {
40954105
if (p.rows()<3) {
@@ -4099,22 +4109,32 @@ namespace BasisClasses
40994109
throw std::runtime_error(msg.str());
41004110
}
41014111

4102-
for (int n=0; n<p.cols(); n++) {
4112+
if (p.rows() == 6) haveVel = true;
41034113

4114+
for (int n=0; n<p.cols(); n++) {
4115+
41044116
if (n % numprocs==myid or not RoundRobin) {
41054117

4118+
for (int k=0; k<3; k++) {
4119+
pp(k) = p(k, n);
4120+
if (haveVel) vv(k) = p(k+3, n);
4121+
}
4122+
4123+
pp = coefrot * (pp - coefctr);
4124+
41064125
bool use = true;
4126+
41074127
if (ftor) {
4108-
for (int k=0; k<3; k++) p1[k] = p(k, n);
4128+
if (haveVel) vv = coefrot * vv;
41094129
use = ftor(m(n), p1, v1, coefindx);
41104130
} else {
41114131
use = true;
41124132
}
41134133
coefindx++;
41144134

4115-
if (use) accumulate(p(0, n)-coefctr[0],
4116-
p(1, n)-coefctr[1],
4117-
p(2, n)-coefctr[2], m(n));
4135+
if (use) {
4136+
accumulate(pp(0), pp(1), pp(2), m(n));
4137+
}
41184138
}
41194139
}
41204140

@@ -4127,22 +4147,32 @@ namespace BasisClasses
41274147
throw std::runtime_error(msg.str());
41284148
}
41294149

4150+
if (p.cols() == 6) haveVel = true;
4151+
4152+
41304153
for (int n=0; n<p.rows(); n++) {
41314154

41324155
if (n % numprocs==myid or not RoundRobin) {
41334156

4157+
for (int k=0; k<3; k++) {
4158+
pp(k) = p(n, k);
4159+
if (haveVel) vv(k) = p(n, k+3);
4160+
}
4161+
4162+
pp = coefrot * (pp - coefctr);
4163+
41344164
bool use = true;
41354165
if (ftor) {
4136-
for (int k=0; k<3; k++) p1[k] = p(n, k);
4166+
if (haveVel) vv = coefrot * vv;
41374167
use = ftor(m(n), p1, v1, coefindx);
41384168
} else {
41394169
use = true;
41404170
}
41414171
coefindx++;
41424172

4143-
if (use) accumulate(p(n, 0)-coefctr[0],
4144-
p(n, 1)-coefctr[1],
4145-
p(n, 2)-coefctr[2], m(n));
4173+
if (use) {
4174+
accumulate(pp(0), pp(1), pp(2), m(n));
4175+
}
41464176
}
41474177
}
41484178
}
@@ -4159,10 +4189,10 @@ namespace BasisClasses
41594189
// Generate coefficients from a phase-space table
41604190
//
41614191
CoefClasses::CoefStrPtr BiorthBasis::createFromArray
4162-
(Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector<double> ctr,
4163-
bool RoundRobin, bool PosVelRows)
4192+
(Eigen::VectorXd& m, RowMatrixXd& p, double time, Eigen::Vector3d ctr,
4193+
RowMatrix3d rot, bool RoundRobin, bool PosVelRows)
41644194
{
4165-
initFromArray(ctr);
4195+
initFromArray(ctr, rot);
41664196
addFromArray(m, p, RoundRobin, PosVelRows);
41674197
return makeFromArray(time);
41684198
}
@@ -4180,13 +4210,20 @@ namespace BasisClasses
41804210
auto ctr = basis->getCenter();
41814211
if (basis->usingNonInertial()) ctr = {0, 0, 0};
41824212

4213+
// Get rotation matrix
4214+
//
4215+
auto rot = basis->getRotation();
4216+
41834217
// Get fields
41844218
//
41854219
int rows = accel.rows();
41864220
for (int n=0; n<rows; n++) {
4187-
auto v = basis->getFields(ps(n, 0) - ctr[0],
4188-
ps(n, 1) - ctr[1],
4189-
ps(n, 2) - ctr[2]);
4221+
Eigen::Vector3d pp;
4222+
for (int k=0; k<3; k++) pp(k) = ps(n, k) - ctr(k);
4223+
pp = rot * pp;
4224+
4225+
auto v = basis->getFields(pp(0), pp(1), pp(2));
4226+
41904227
// First 6 fields are density and potential, followed by acceleration
41914228
for (int k=0; k<3; k++) accel(n, k) += v[6+k] - basis->pseudo(k);
41924229
}
@@ -4271,11 +4308,18 @@ namespace BasisClasses
42714308

42724309
// Interpolate center
42734310
//
4274-
if (coefsA->ctr.size() and coefsB->ctr.size()) {
4275-
newcoef->ctr.resize(3);
4276-
for (int k=0; k<3; k++)
4277-
newcoef->ctr[k] = a * coefsA->ctr[k] + b * coefsB->ctr[k];
4278-
}
4311+
newcoef->ctr = a * coefsA->ctr + b * coefsB->ctr;
4312+
4313+
// Interpolate rotation matrix followed by unitarization
4314+
//
4315+
RowMatrix3d newrot = a * coefsA->rot + b * coefsB->rot;
4316+
4317+
// Closest unitary matrix in the Frobenius norm sense
4318+
//
4319+
Eigen::BDCSVD<RowMatrix3d> svd
4320+
(newrot, Eigen::ComputeFullU | Eigen::ComputeFullV);
4321+
4322+
newcoef->rot = svd.matrixU() * svd.matrixV().adjoint();
42794323

42804324
// Install coefficients
42814325
//

0 commit comments

Comments
 (0)