Skip to content

Commit c57f653

Browse files
author
Martin D. Weinberg
committed
Preliminary implementation of rotation matrix enhancement
1 parent cb6c7c1 commit c57f653

11 files changed

Lines changed: 208 additions & 43 deletions

‎expui/BasisFactory.H‎

Lines changed: 14 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,9 @@ namespace BasisClasses
6868
//! The expansion center
6969
std::vector<double> coefctr;
7070

71+
//! Rotation matrix
72+
Eigen::Matrix3d coefrot;
73+
7174
//! Contains contructor and BFE parameter database
7275
YAML::Node node, conf;
7376

@@ -151,8 +154,10 @@ namespace BasisClasses
151154
virtual ~Basis(void) {}
152155

153156
//! Set the expansion center
154-
void setCenter(std::vector<double> center)
155-
{ coefctr = center; }
157+
void setCenter(std::vector<double> center) { coefctr = center; }
158+
159+
//! Set the rotation matrix
160+
void setRotation(Eigen::Matrix3d rot) { coefrot = rot; }
156161

157162
//! Evaluate basis in desired coordinates
158163
virtual std::vector<double>
@@ -213,11 +218,12 @@ namespace BasisClasses
213218

214219
//! Generate coeffients from a particle reader
215220
virtual CoefClasses::CoefStrPtr
216-
createFromReader(PR::PRptr reader, std::vector<double> ctr) = 0;
221+
createFromReader(PR::PRptr reader, std::vector<double> ctr,
222+
Eigen::Matrix3d rot) = 0;
217223

218224
//! Generate coefficients from a phase-space table
219225
virtual void
220-
initFromArray(std::vector<double> ctr) = 0;
226+
initFromArray(std::vector<double> ctr, Eigen::Matrix3d rot) = 0;
221227

222228
//! Accumulate coefficient contributions from arrays
223229
virtual void
@@ -229,6 +235,7 @@ namespace BasisClasses
229235
CoefClasses::CoefStrPtr createFromArray
230236
(Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0,
231237
std::vector<double> center={0.0, 0.0, 0.0},
238+
Eigen::Matrix3d rot=Eigen::Matrix3d::Identity(),
232239
bool roundrobin=true, bool posvelrows=false);
233240

234241
//! Create and the coefficients from the array accumulation with the
@@ -276,6 +283,9 @@ namespace BasisClasses
276283

277284
//! Get the basis expansion center
278285
std::vector<double> getCenter() { return coefctr; }
286+
287+
//! Get the basis expansion center
288+
Eigen::Matrix3d getRotation() { return coefrot; }
279289
};
280290

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

‎expui/BasisFactory.cc‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -276,9 +276,9 @@ namespace BasisClasses
276276
//
277277
CoefClasses::CoefStrPtr Basis::createFromArray
278278
(Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector<double> ctr,
279-
bool roundrobin, bool posvelrows)
279+
Eigen::Matrix3d 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: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -94,13 +94,15 @@ namespace BasisClasses
9494
//! Generate coeffients from a particle reader and optional center
9595
//! location for the expansion
9696
CoefClasses::CoefStrPtr createFromReader
97-
(PR::PRptr reader, std::vector<double> center={0.0, 0.0, 0.0});
97+
(PR::PRptr reader, std::vector<double> center={0.0, 0.0, 0.0},
98+
Eigen::Matrix3d rot=Eigen::Matrix3d::Identity());
9899

99100
//! Generate coeffients from an array and optional center location
100101
//! for the expansion
101102
CoefClasses::CoefStrPtr createFromArray
102103
(Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0,
103104
std::vector<double> center={0.0, 0.0, 0.0},
105+
Eigen::Matrix3d rot=Eigen::Matrix3d::Identity(),
104106
bool roundrobin=true, bool posvelrows=false);
105107

106108
//! Generate coeffients from an array and optional center location
@@ -110,7 +112,8 @@ namespace BasisClasses
110112
//! Initialize accumulating coefficients from arrays with an optional
111113
//! center vector. This is called once to initialize the accumulation.
112114
void initFromArray
113-
(std::vector<double> center={0.0, 0.0, 0.0});
115+
(std::vector<double> center={0.0, 0.0, 0.0},
116+
Eigen::Matrix3d rot=Eigen::Matrix3d::Identity());
114117

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

‎expui/BiorthBasis.cc‎

Lines changed: 59 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -3871,7 +3871,7 @@ namespace BasisClasses
38713871

38723872
// Generate coeffients from a particle reader
38733873
CoefClasses::CoefStrPtr BiorthBasis::createFromReader
3874-
(PR::PRptr reader, std::vector<double> ctr)
3874+
(PR::PRptr reader, std::vector<double> ctr, Eigen::Matrix3d rot)
38753875
{
38763876
CoefClasses::CoefStrPtr coef;
38773877

@@ -3901,6 +3901,10 @@ namespace BasisClasses
39013901
//
39023902
if (addCenter) coef->ctr = ctr;
39033903

3904+
// Add the rotation matrix
3905+
//
3906+
coef->rot = rot;
3907+
39043908
std::vector<double> pp(3), vv(3);
39053909

39063910
reset_coefs();
@@ -3916,18 +3920,21 @@ namespace BasisClasses
39163920
use = true;
39173921
}
39183922

3919-
if (use) accumulate(p->pos[0]-ctr[0],
3920-
p->pos[1]-ctr[1],
3921-
p->pos[2]-ctr[2],
3922-
p->mass);
3923+
if (use) {
3924+
Eigen::Vector3d pp;
3925+
for (int k=0; k<3; k++) pp(k) = p->pos[k] - coefctr[k];
3926+
pp = coefrot * pp;
3927+
3928+
accumulate(pp(0), pp(1), pp(2), p->mass);
3929+
}
39233930
}
39243931
make_coefs();
39253932
load_coefs(coef, reader->CurrentTime());
39263933
return coef;
39273934
}
39283935

39293936
// Generate coefficients from a phase-space table
3930-
void BiorthBasis::initFromArray(std::vector<double> ctr)
3937+
void BiorthBasis::initFromArray(std::vector<double> ctr, Eigen::Matrix3d rot)
39313938
{
39323939
if (name.compare("sphereSL") == 0)
39333940
coefret = std::make_shared<CoefClasses::SphStruct>();
@@ -3957,6 +3964,9 @@ namespace BasisClasses
39573964
//
39583965
coefctr = ctr;
39593966

3967+
// Register the rotation matrix
3968+
coefret->rot = rot;
3969+
39603970
// Clean up for accumulation
39613971
//
39623972
reset_coefs();
@@ -4027,9 +4037,13 @@ namespace BasisClasses
40274037
}
40284038
coefindx++;
40294039

4030-
if (use) accumulate(p(0, n)-coefctr[0],
4031-
p(1, n)-coefctr[1],
4032-
p(2, n)-coefctr[2], m(n));
4040+
if (use) {
4041+
Eigen::Vector3d pp;
4042+
for (int k=0; k<3; k++) pp(k) = p(k, n) - coefctr[k];
4043+
pp = coefrot * pp;
4044+
4045+
accumulate(pp(0), pp(1), pp(2), m(n));
4046+
}
40334047
}
40344048
}
40354049

@@ -4055,9 +4069,13 @@ namespace BasisClasses
40554069
}
40564070
coefindx++;
40574071

4058-
if (use) accumulate(p(n, 0)-coefctr[0],
4059-
p(n, 1)-coefctr[1],
4060-
p(n, 2)-coefctr[2], m(n));
4072+
if (use) {
4073+
Eigen::Vector3d pp;
4074+
for (int k=0; k<3; k++) pp(k) = p(k, n) - coefctr[0];
4075+
pp = coefrot * pp;
4076+
4077+
accumulate(pp(0), pp(1), pp(2), m(n));
4078+
}
40614079
}
40624080
}
40634081
}
@@ -4075,9 +4093,9 @@ namespace BasisClasses
40754093
//
40764094
CoefClasses::CoefStrPtr BiorthBasis::createFromArray
40774095
(Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector<double> ctr,
4078-
bool RoundRobin, bool PosVelRows)
4096+
Eigen::Matrix3d rot, bool RoundRobin, bool PosVelRows)
40794097
{
4080-
initFromArray(ctr);
4098+
initFromArray(ctr, rot);
40814099
addFromArray(m, p, RoundRobin, PosVelRows);
40824100
return makeFromArray(time);
40834101
}
@@ -4095,13 +4113,20 @@ namespace BasisClasses
40954113
auto ctr = basis->getCenter();
40964114
if (basis->usingNonInertial()) ctr = {0, 0, 0};
40974115

4116+
// Get rotation matrix
4117+
//
4118+
auto rot = basis->getRotation();
4119+
40984120
// Get fields
40994121
//
41004122
int rows = accel.rows();
41014123
for (int n=0; n<rows; n++) {
4102-
auto v = basis->getFields(ps(n, 0) - ctr[0],
4103-
ps(n, 1) - ctr[1],
4104-
ps(n, 2) - ctr[2]);
4124+
Eigen::Vector3d pp;
4125+
for (int k=0; k<3; k++) pp(k) = ps(n, k) - ctr[k];
4126+
pp = rot * pp;
4127+
4128+
auto v = basis->getFields(pp(0), pp(1), pp(2));
4129+
41054130
// First 6 fields are density and potential, followed by acceleration
41064131
for (int k=0; k<3; k++) accel(n, k) += v[6+k] - basis->pseudo(k);
41074132
}
@@ -4192,6 +4217,23 @@ namespace BasisClasses
41924217
newcoef->ctr[k] = a * coefsA->ctr[k] + b * coefsB->ctr[k];
41934218
}
41944219

4220+
// Interpolate rotation matrix followed by unitarization
4221+
//
4222+
Eigen::Matrix3d newrot;
4223+
for (int k=0; k<9; k++) {
4224+
newrot.data()[k] =
4225+
a * coefsA->rot.data()[k] + b * coefsB->rot.data()[k];
4226+
}
4227+
4228+
// Closest unitary matrix in the Frobenius norm sense
4229+
//
4230+
Eigen::BDCSVD<Eigen::Matrix3d> svd(newrot,
4231+
Eigen::ComputeFullU | Eigen::ComputeFullV);
4232+
auto U = svd.matrixU();
4233+
auto V = svd.matrixV();
4234+
4235+
newcoef->rot = U * V.adjoint();
4236+
41954237
// Install coefficients
41964238
//
41974239
basis->set_coefs(newcoef);

‎expui/CoefStruct.H‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,9 @@ namespace CoefClasses
2929
//! Center data
3030
std::vector<double> ctr= {0.0, 0.0, 0.0};
3131

32+
//! Orientation data
33+
Eigen::Matrix3d rot= Eigen::Matrix3d::Identity();
34+
3235
//! Destructor
3336
virtual ~CoefStruct() {}
3437

@@ -73,6 +76,12 @@ namespace CoefClasses
7376
//! Read-only access to center (no copy)
7477
std::vector<double> getCenter() { return ctr; }
7578

79+
//! Set new rotation matrix
80+
void setRotation(Eigen::Matrix3d& ROT) { rot = ROT; }
81+
82+
//! Read-only access to orientation
83+
Eigen::Matrix3d getRotation() { return rot; }
84+
7685
//! Set coefficient time (no copy)
7786
void setTime(double& STORE)
7887
{

‎expui/Coefficients.cc‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -118,6 +118,13 @@ namespace CoefClasses
118118
stanza.getAttribute("Center").read(ctr);
119119
}
120120

121+
// Check for rotation matrix
122+
//
123+
Eigen::Matrix3d rot = Eigen::Matrix3d::Identity();
124+
if (stanza.hasAttribute("Rotation")) {
125+
stanza.getAttribute("Rotation").read(rot);
126+
}
127+
121128
if (Time < Tmin or Time > Tmax) continue;
122129

123130
auto in = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
@@ -143,6 +150,7 @@ namespace CoefClasses
143150

144151
if (ctr.size()) coef->ctr = ctr;
145152

153+
coef->rot = rot;
146154
coef->lmax = Lmax;
147155
coef->nmax = Nmax;
148156
coef->time = Time;
@@ -318,6 +326,13 @@ namespace CoefClasses
318326
stanza.getAttribute("Center").read(ctr);
319327
}
320328

329+
// Check for rotation matrix
330+
//
331+
Eigen::Matrix3d rot = Eigen::Matrix3d::Identity();
332+
if (stanza.hasAttribute("Rotation")) {
333+
stanza.getAttribute("Rotation").read(rot);
334+
}
335+
321336
if (Time < Tmin or Time > Tmax) continue;
322337

323338
std::array<long int, 3> shape;
@@ -331,6 +346,7 @@ namespace CoefClasses
331346

332347
if (ctr.size()) coef->ctr = ctr;
333348

349+
coef->rot = rot;
334350
coef->nfld = Nfld;
335351
coef->lmax = Lmax;
336352
coef->nmax = Nmax;
@@ -409,6 +425,14 @@ namespace CoefClasses
409425
stanza.getAttribute("Center").read(ctr);
410426
}
411427

428+
// Check for rotation matrix
429+
//
430+
Eigen::Matrix3d rot = Eigen::Matrix3d::Identity();
431+
if (stanza.hasAttribute("Rotation")) {
432+
stanza.getAttribute("Rotation").read(rot);
433+
}
434+
435+
412436
if (Time < Tmin or Time > Tmax) continue;
413437

414438
std::array<long int, 3> shape;
@@ -422,6 +446,7 @@ namespace CoefClasses
422446

423447
if (ctr.size()) coef->ctr = ctr;
424448

449+
coef->rot = rot;
425450
coef->nfld = Nfld;
426451
coef->mmax = Mmax;
427452
coef->nmax = Nmax;
@@ -665,6 +690,11 @@ namespace CoefClasses
665690
if (C->ctr.size()>0)
666691
stanza.createAttribute<double>("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr);
667692

693+
// Add a rotation matrix attribute
694+
//
695+
Eigen::Matrix3d rot = C->getRotation();
696+
stanza.createAttribute<double>("Rotation", HighFive::DataSpace::From(rot)).write(rot);
697+
668698
// Index counters
669699
//
670700
unsigned I = 0, L = 0;
@@ -1056,6 +1086,11 @@ namespace CoefClasses
10561086
if (C->ctr.size()>0)
10571087
stanza.createAttribute<double>("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr);
10581088

1089+
// Add a rotation matrix attribute
1090+
//
1091+
Eigen::Matrix3d rot = C->getRotation();
1092+
stanza.createAttribute<Eigen::Matrix3d>("Rotation", HighFive::DataSpace::From(rot)).write(rot);
1093+
10591094
// Add coefficient data
10601095
//
10611096
Eigen::MatrixXcd out(*C->coefs);
@@ -2921,6 +2956,11 @@ namespace CoefClasses
29212956
stanza.createAttribute<double>("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr);
29222957

29232958

2959+
// Add a rotation matrix attribute
2960+
//
2961+
Eigen::Matrix3d rot = C->getRotation();
2962+
stanza.createAttribute<Eigen::Matrix3d>("Rotation", HighFive::DataSpace::From(rot)).write(rot);
2963+
29242964
// Coefficient size (allow Eigen::Tensor to be easily recontructed from metadata)
29252965
//
29262966
const auto& d = C->coefs->dimensions();

‎expui/FieldBasis.H‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -135,7 +135,8 @@ namespace BasisClasses
135135
//! Generate coeffients from a particle reader and optional center
136136
//! location for the expansion
137137
CoefClasses::CoefStrPtr createFromReader
138-
(PR::PRptr reader, std::vector<double> center={0.0, 0.0, 0.0});
138+
(PR::PRptr reader, std::vector<double> center={0.0, 0.0, 0.0},
139+
Eigen::Matrix3d rot=Eigen::Matrix3d::Identity());
139140

140141

141142
//! Generate coeffients from an array and optional center location
@@ -145,7 +146,9 @@ namespace BasisClasses
145146
//! Initialize accumulating coefficients from arrays with an optional
146147
//! center vector. This is called once to initialize the accumulation.
147148
void initFromArray
148-
(std::vector<double> center={0.0, 0.0, 0.0});
149+
(std::vector<double> center={0.0, 0.0, 0.0},
150+
Eigen::Matrix3d rot=Eigen::Matrix3d::Identity());
151+
149152

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

0 commit comments

Comments
 (0)