From f6f2f3f28c71982d18819754fecd68a8e30bebd5 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 18 Jul 2025 12:46:41 -0400 Subject: [PATCH 01/18] Add inner/outer boundary restrictions for SL solution to prevent underflow while computing basis elements with large harmonic orders --- exputil/SLGridMP2.cc | 102 +++++++++++++++++++++++++++++++++++++------ include/SLGridMP2.H | 12 +++++ 2 files changed, 100 insertions(+), 14 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 23f485f42..bb358fd0b 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -120,6 +120,13 @@ double sphdens(double r) return 4.0*M_PI * model->get_density(r); } +// Use entire grid for for l=Lswitch: rmin=rmap/rfac, rmax=rmap*rfac with rfac=pow(10, +// Lalpha/l) +double SLGridSph::Lalpha = 100.0; + void SLGridSph::bomb(string oops) { @@ -1109,10 +1116,25 @@ void SLGridSph::compute_table(struct TableSph* table, int l) if (myid==0) VERBOSE = SLEDGE_VERBOSE; #endif - cons[6] = rmin; - cons[7] = rmax; + // Inner radial boundary heuristic + // + int Nlo = 0, Nhi = numr; + if (l>Lswitch) { + double Rfac = std::pow(10.0, Lalpha/l); + double Rmin = std::max(rmin, rmap/Rfac); + double Rmax = std::min(rmax, rmap*Rfac); + + // Find index boundaries in r grid + auto lower = std::lower_bound(r.data(), r.data() + r.size(), Rmin); + Nlo = std::distance(r.data(), lower); + auto upper = std::upper_bound(r.data(), r.data() + r.size(), Rmax); + Nhi = std::distance(r.data(), upper); + } + + cons[6] = r[Nlo]; + cons[7] = r[Nhi-1]; L2 = l*(l+1); - NUM = numr; + NUM = Nhi - Nlo; N = nmax; // integer iflag[nmax], invec[nmax+3]; @@ -1164,7 +1186,7 @@ void SLGridSph::compute_table(struct TableSph* table, int l) // // Output mesh // - for (int i=0; i -10) { + std::cout << "# order=" << i + << " error: " << sledge_error(iflag[i]) + << std::endl << "#" << std::endl; + + std::cout << std::setw(14) << "x" + << std::setw(25) << "u(x)" + << std::setw(25) << "(pu`)(x)" + << std::endl; + int k = NUM*i; + for (int j=0; jev[i] = ev[i]; table->ef.resize(N, numr); + table->ef.setZero(); // Choose sign conventions for the ef table // @@ -1282,9 +1332,11 @@ void SLGridSph::compute_table(struct TableSph* table, int l) if (ef[j*NUM+nfid]<0.0) sgn(j) = -1; } - for (int i=0; ief(j, i) = ef[j*NUM+i] * sgn(j); + table->ef(j, i+Nlo) = ef[j*NUM+i] * sgn(j); } table->l = l; @@ -1374,11 +1426,26 @@ void SLGridSph::compute_table_worker(void) if (tbdbg) std::cout << "Worker " << mpi_myid << ": ordered to compute l = " << L << "" << std::endl; + // Inner radial boundary heuristic + // + int Nlo = 0, Nhi = numr; + if (L>Lswitch) { + double Rfac = std::pow(10.0, Lalpha/L); + double Rmin = std::max(rmin, rmap/Rfac); + double Rmax = std::min(rmax, rmap*Rfac); + + // Find index boundaries in r grid + auto lower = std::lower_bound(r.data(), r.data() + r.size(), Rmin); + Nlo = std::distance(r.data(), lower); + auto upper = std::upper_bound(r.data(), r.data() + r.size(), Rmax); + Nhi = std::distance(r.data(), upper); + } + cons[0] = cons[1] = cons[2] = cons[3] = cons[4] = cons[5] = 0.0; - cons[6] = rmin; - cons[7] = rmax; + cons[6] = r[Nlo]; + cons[7] = r[Nhi-1]; L2 = L*(L+1); - NUM = numr; + NUM = Nhi - Nlo; N = nmax; // integer iflag[nmax], invec[nmax+3]; @@ -1422,7 +1489,7 @@ void SLGridSph::compute_table_worker(void) // // Output mesh // - for (int i=0; i -10) { cout << std::setw(14) << "x" diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index 5bb6888fe..052dfb6ab 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -94,6 +94,18 @@ private: //! Cache versioning inline static const std::string Version = "1.0"; + //@{ + /** Rmin heuristic for grid construction to prevent underflow issues + in Sturm-Liouville solver: sledge */ + + //! Use entire grid for l=Lswitch: + //! rmin=rmap/rfac, rmax=rmap*rfac with rfac=pow(10, Lalpha/l) + static double Lalpha; + //@} + public: //! Flag for MPI enabled (default: 0=off) From c57f6539e6eed6c0bfff132404ea27a6bc78fc23 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 29 Jul 2025 13:36:43 -0400 Subject: [PATCH 02/18] Preliminary implementation of rotation matrix enhancement --- expui/BasisFactory.H | 18 +++++++--- expui/BasisFactory.cc | 4 +-- expui/BiorthBasis.H | 7 ++-- expui/BiorthBasis.cc | 76 ++++++++++++++++++++++++++++++++---------- expui/CoefStruct.H | 9 +++++ expui/Coefficients.cc | 40 ++++++++++++++++++++++ expui/FieldBasis.H | 7 ++-- expui/FieldBasis.cc | 4 +-- pyEXP/BasisWrappers.cc | 36 ++++++++++++-------- pyEXP/CoefWrappers.cc | 40 ++++++++++++++++++++++ src/Component.H | 10 +++++- 11 files changed, 208 insertions(+), 43 deletions(-) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index 1598f3f7d..c16fffab5 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -68,6 +68,9 @@ namespace BasisClasses //! The expansion center std::vector coefctr; + //! Rotation matrix + Eigen::Matrix3d coefrot; + //! Contains contructor and BFE parameter database YAML::Node node, conf; @@ -151,8 +154,10 @@ namespace BasisClasses virtual ~Basis(void) {} //! Set the expansion center - void setCenter(std::vector center) - { coefctr = center; } + void setCenter(std::vector center) { coefctr = center; } + + //! Set the rotation matrix + void setRotation(Eigen::Matrix3d rot) { coefrot = rot; } //! Evaluate basis in desired coordinates virtual std::vector @@ -213,11 +218,12 @@ namespace BasisClasses //! Generate coeffients from a particle reader virtual CoefClasses::CoefStrPtr - createFromReader(PR::PRptr reader, std::vector ctr) = 0; + createFromReader(PR::PRptr reader, std::vector ctr, + Eigen::Matrix3d rot) = 0; //! Generate coefficients from a phase-space table virtual void - initFromArray(std::vector ctr) = 0; + initFromArray(std::vector ctr, Eigen::Matrix3d rot) = 0; //! Accumulate coefficient contributions from arrays virtual void @@ -229,6 +235,7 @@ namespace BasisClasses CoefClasses::CoefStrPtr createFromArray (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0, std::vector center={0.0, 0.0, 0.0}, + Eigen::Matrix3d rot=Eigen::Matrix3d::Identity(), bool roundrobin=true, bool posvelrows=false); //! Create and the coefficients from the array accumulation with the @@ -276,6 +283,9 @@ namespace BasisClasses //! Get the basis expansion center std::vector getCenter() { return coefctr; } + + //! Get the basis expansion center + Eigen::Matrix3d getRotation() { return coefrot; } }; using BasisPtr = std::shared_ptr; diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index b410c7861..0e486933b 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -276,9 +276,9 @@ namespace BasisClasses // CoefClasses::CoefStrPtr Basis::createFromArray (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector ctr, - bool roundrobin, bool posvelrows) + Eigen::Matrix3d rot, bool roundrobin, bool posvelrows) { - initFromArray(ctr); + initFromArray(ctr, rot); addFromArray(m, p, roundrobin, posvelrows); return makeFromArray(time); } diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index b8e1609bf..7e2f70d12 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -94,13 +94,15 @@ namespace BasisClasses //! Generate coeffients from a particle reader and optional center //! location for the expansion CoefClasses::CoefStrPtr createFromReader - (PR::PRptr reader, std::vector center={0.0, 0.0, 0.0}); + (PR::PRptr reader, std::vector center={0.0, 0.0, 0.0}, + Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); //! Generate coeffients from an array and optional center location //! for the expansion CoefClasses::CoefStrPtr createFromArray (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0, std::vector center={0.0, 0.0, 0.0}, + Eigen::Matrix3d rot=Eigen::Matrix3d::Identity(), bool roundrobin=true, bool posvelrows=false); //! Generate coeffients from an array and optional center location @@ -110,7 +112,8 @@ namespace BasisClasses //! Initialize accumulating coefficients from arrays with an optional //! center vector. This is called once to initialize the accumulation. void initFromArray - (std::vector center={0.0, 0.0, 0.0}); + (std::vector center={0.0, 0.0, 0.0}, + Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); //! Initialize accumulating coefficients from arrays. This is //! called once to initialize the accumulation. diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 203f5f115..5f7a061a4 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3871,7 +3871,7 @@ namespace BasisClasses // Generate coeffients from a particle reader CoefClasses::CoefStrPtr BiorthBasis::createFromReader - (PR::PRptr reader, std::vector ctr) + (PR::PRptr reader, std::vector ctr, Eigen::Matrix3d rot) { CoefClasses::CoefStrPtr coef; @@ -3901,6 +3901,10 @@ namespace BasisClasses // if (addCenter) coef->ctr = ctr; + // Add the rotation matrix + // + coef->rot = rot; + std::vector pp(3), vv(3); reset_coefs(); @@ -3916,10 +3920,13 @@ namespace BasisClasses use = true; } - if (use) accumulate(p->pos[0]-ctr[0], - p->pos[1]-ctr[1], - p->pos[2]-ctr[2], - p->mass); + if (use) { + Eigen::Vector3d pp; + for (int k=0; k<3; k++) pp(k) = p->pos[k] - coefctr[k]; + pp = coefrot * pp; + + accumulate(pp(0), pp(1), pp(2), p->mass); + } } make_coefs(); load_coefs(coef, reader->CurrentTime()); @@ -3927,7 +3934,7 @@ namespace BasisClasses } // Generate coefficients from a phase-space table - void BiorthBasis::initFromArray(std::vector ctr) + void BiorthBasis::initFromArray(std::vector ctr, Eigen::Matrix3d rot) { if (name.compare("sphereSL") == 0) coefret = std::make_shared(); @@ -3957,6 +3964,9 @@ namespace BasisClasses // coefctr = ctr; + // Register the rotation matrix + coefret->rot = rot; + // Clean up for accumulation // reset_coefs(); @@ -4027,9 +4037,13 @@ namespace BasisClasses } coefindx++; - if (use) accumulate(p(0, n)-coefctr[0], - p(1, n)-coefctr[1], - p(2, n)-coefctr[2], m(n)); + if (use) { + Eigen::Vector3d pp; + for (int k=0; k<3; k++) pp(k) = p(k, n) - coefctr[k]; + pp = coefrot * pp; + + accumulate(pp(0), pp(1), pp(2), m(n)); + } } } @@ -4055,9 +4069,13 @@ namespace BasisClasses } coefindx++; - if (use) accumulate(p(n, 0)-coefctr[0], - p(n, 1)-coefctr[1], - p(n, 2)-coefctr[2], m(n)); + if (use) { + Eigen::Vector3d pp; + for (int k=0; k<3; k++) pp(k) = p(k, n) - coefctr[0]; + pp = coefrot * pp; + + accumulate(pp(0), pp(1), pp(2), m(n)); + } } } } @@ -4075,9 +4093,9 @@ namespace BasisClasses // CoefClasses::CoefStrPtr BiorthBasis::createFromArray (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector ctr, - bool RoundRobin, bool PosVelRows) + Eigen::Matrix3d rot, bool RoundRobin, bool PosVelRows) { - initFromArray(ctr); + initFromArray(ctr, rot); addFromArray(m, p, RoundRobin, PosVelRows); return makeFromArray(time); } @@ -4095,13 +4113,20 @@ namespace BasisClasses auto ctr = basis->getCenter(); if (basis->usingNonInertial()) ctr = {0, 0, 0}; + // Get rotation matrix + // + auto rot = basis->getRotation(); + // Get fields // int rows = accel.rows(); for (int n=0; ngetFields(ps(n, 0) - ctr[0], - ps(n, 1) - ctr[1], - ps(n, 2) - ctr[2]); + Eigen::Vector3d pp; + for (int k=0; k<3; k++) pp(k) = ps(n, k) - ctr[k]; + pp = rot * pp; + + auto v = basis->getFields(pp(0), pp(1), pp(2)); + // First 6 fields are density and potential, followed by acceleration for (int k=0; k<3; k++) accel(n, k) += v[6+k] - basis->pseudo(k); } @@ -4192,6 +4217,23 @@ namespace BasisClasses newcoef->ctr[k] = a * coefsA->ctr[k] + b * coefsB->ctr[k]; } + // Interpolate rotation matrix followed by unitarization + // + Eigen::Matrix3d newrot; + for (int k=0; k<9; k++) { + newrot.data()[k] = + a * coefsA->rot.data()[k] + b * coefsB->rot.data()[k]; + } + + // Closest unitary matrix in the Frobenius norm sense + // + Eigen::BDCSVD svd(newrot, + Eigen::ComputeFullU | Eigen::ComputeFullV); + auto U = svd.matrixU(); + auto V = svd.matrixV(); + + newcoef->rot = U * V.adjoint(); + // Install coefficients // basis->set_coefs(newcoef); diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index 26d011e28..052a9e4e0 100644 --- a/expui/CoefStruct.H +++ b/expui/CoefStruct.H @@ -29,6 +29,9 @@ namespace CoefClasses //! Center data std::vector ctr= {0.0, 0.0, 0.0}; + //! Orientation data + Eigen::Matrix3d rot= Eigen::Matrix3d::Identity(); + //! Destructor virtual ~CoefStruct() {} @@ -73,6 +76,12 @@ namespace CoefClasses //! Read-only access to center (no copy) std::vector getCenter() { return ctr; } + //! Set new rotation matrix + void setRotation(Eigen::Matrix3d& ROT) { rot = ROT; } + + //! Read-only access to orientation + Eigen::Matrix3d getRotation() { return rot; } + //! Set coefficient time (no copy) void setTime(double& STORE) { diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index b01cb0c04..10e2fe233 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -118,6 +118,13 @@ namespace CoefClasses stanza.getAttribute("Center").read(ctr); } + // Check for rotation matrix + // + Eigen::Matrix3d rot = Eigen::Matrix3d::Identity(); + if (stanza.hasAttribute("Rotation")) { + stanza.getAttribute("Rotation").read(rot); + } + if (Time < Tmin or Time > Tmax) continue; auto in = stanza.getDataSet("coefficients").read(); @@ -143,6 +150,7 @@ namespace CoefClasses if (ctr.size()) coef->ctr = ctr; + coef->rot = rot; coef->lmax = Lmax; coef->nmax = Nmax; coef->time = Time; @@ -318,6 +326,13 @@ namespace CoefClasses stanza.getAttribute("Center").read(ctr); } + // Check for rotation matrix + // + Eigen::Matrix3d rot = Eigen::Matrix3d::Identity(); + if (stanza.hasAttribute("Rotation")) { + stanza.getAttribute("Rotation").read(rot); + } + if (Time < Tmin or Time > Tmax) continue; std::array shape; @@ -331,6 +346,7 @@ namespace CoefClasses if (ctr.size()) coef->ctr = ctr; + coef->rot = rot; coef->nfld = Nfld; coef->lmax = Lmax; coef->nmax = Nmax; @@ -409,6 +425,14 @@ namespace CoefClasses stanza.getAttribute("Center").read(ctr); } + // Check for rotation matrix + // + Eigen::Matrix3d rot = Eigen::Matrix3d::Identity(); + if (stanza.hasAttribute("Rotation")) { + stanza.getAttribute("Rotation").read(rot); + } + + if (Time < Tmin or Time > Tmax) continue; std::array shape; @@ -422,6 +446,7 @@ namespace CoefClasses if (ctr.size()) coef->ctr = ctr; + coef->rot = rot; coef->nfld = Nfld; coef->mmax = Mmax; coef->nmax = Nmax; @@ -665,6 +690,11 @@ namespace CoefClasses if (C->ctr.size()>0) stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); + // Add a rotation matrix attribute + // + Eigen::Matrix3d rot = C->getRotation(); + stanza.createAttribute("Rotation", HighFive::DataSpace::From(rot)).write(rot); + // Index counters // unsigned I = 0, L = 0; @@ -1056,6 +1086,11 @@ namespace CoefClasses if (C->ctr.size()>0) stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); + // Add a rotation matrix attribute + // + Eigen::Matrix3d rot = C->getRotation(); + stanza.createAttribute("Rotation", HighFive::DataSpace::From(rot)).write(rot); + // Add coefficient data // Eigen::MatrixXcd out(*C->coefs); @@ -2921,6 +2956,11 @@ namespace CoefClasses stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); + // Add a rotation matrix attribute + // + Eigen::Matrix3d rot = C->getRotation(); + stanza.createAttribute("Rotation", HighFive::DataSpace::From(rot)).write(rot); + // Coefficient size (allow Eigen::Tensor to be easily recontructed from metadata) // const auto& d = C->coefs->dimensions(); diff --git a/expui/FieldBasis.H b/expui/FieldBasis.H index f22893fe2..432706c19 100644 --- a/expui/FieldBasis.H +++ b/expui/FieldBasis.H @@ -135,7 +135,8 @@ namespace BasisClasses //! Generate coeffients from a particle reader and optional center //! location for the expansion CoefClasses::CoefStrPtr createFromReader - (PR::PRptr reader, std::vector center={0.0, 0.0, 0.0}); + (PR::PRptr reader, std::vector center={0.0, 0.0, 0.0}, + Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); //! Generate coeffients from an array and optional center location @@ -145,7 +146,9 @@ namespace BasisClasses //! Initialize accumulating coefficients from arrays with an optional //! center vector. This is called once to initialize the accumulation. void initFromArray - (std::vector center={0.0, 0.0, 0.0}); + (std::vector center={0.0, 0.0, 0.0}, + Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); + //! Initialize accumulating coefficients from arrays. This is //! called once to initialize the accumulation. diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index ff59c2016..a720afaac 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -539,7 +539,7 @@ namespace BasisClasses // Generate coeffients from a particle reader CoefClasses::CoefStrPtr FieldBasis::createFromReader - (PR::PRptr reader, std::vector ctr) + (PR::PRptr reader, std::vector ctr, Eigen::Matrix3d rot) { CoefClasses::CoefStrPtr coef; @@ -595,7 +595,7 @@ namespace BasisClasses } // Generate coefficients from a phase-space table - void FieldBasis::initFromArray(std::vector ctr) + void FieldBasis::initFromArray(std::vector ctr, Eigen::Matrix3d rot) { if (dof==3) coefret = std::make_shared(); diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 58042c8fb..0f91d0f06 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -292,13 +292,15 @@ void BasisFactoryClasses(py::module &m) } virtual CoefClasses::CoefStrPtr - createFromReader(PR::PRptr reader, std::vector ctr) override { - PYBIND11_OVERRIDE_PURE(CoefClasses::CoefStrPtr, Basis, createFromReader, reader, ctr); + createFromReader(PR::PRptr reader, std::vector ctr, + Eigen::Matrix3d rot) override { + PYBIND11_OVERRIDE_PURE(CoefClasses::CoefStrPtr, Basis, createFromReader, reader, ctr, rot); } - virtual void initFromArray(std::vector ctr) override { - PYBIND11_OVERRIDE_PURE(void, Basis, initFromArray, ctr); + virtual void initFromArray + (std::vector ctr, Eigen::MatrixXd rot) { + PYBIND11_OVERRIDE_PURE(void, Basis, initFromArray, ctr, rot); } virtual void @@ -886,10 +888,10 @@ void BasisFactoryClasses(py::module &m) ) .def("createFromArray", [](BasisClasses::Basis& A, Eigen::VectorXd& mass, RowMatrixXd& ps, - double time, std::vector center, + double time, std::vector center, Eigen::Matrix3d rot, bool roundrobin, bool posvelrows) { - return A.createFromArray(mass, ps, time, center, + return A.createFromArray(mass, ps, time, center, rot, roundrobin, posvelrows); }, R"( @@ -926,6 +928,7 @@ void BasisFactoryClasses(py::module &m) )", py::arg("mass"), py::arg("pos"), py::arg("time"), py::arg("center") = std::vector(3, 0.0), + py::arg("rotation") = Eigen::Matrix3d::Identity(), py::arg("roundrobin") = true, py::arg("posvelrows") = true) .def("makeFromArray", [](BasisClasses::Basis& A, double time) @@ -1127,13 +1130,14 @@ void BasisFactoryClasses(py::module &m) the basis coefficients computed from the particles )", py::arg("reader"), - py::arg("center") = std::vector(3, 0.0)) + py::arg("center") = std::vector(3, 0.0), + py::arg("rotation") = Eigen::Matrix3d::Identity()) .def("createFromArray", [](BasisClasses::BiorthBasis& A, Eigen::VectorXd& mass, RowMatrixXd& pos, - double time, std::vector center, + double time, std::vector center, Eigen::Matrix3d rot, bool roundrobin, bool posvelrows) { - return A.createFromArray(mass, pos, time, center, + return A.createFromArray(mass, pos, time, center, rot, roundrobin, posvelrows); }, R"( @@ -1169,6 +1173,7 @@ void BasisFactoryClasses(py::module &m) )", py::arg("mass"), py::arg("pos"), py::arg("time"), py::arg("center") = std::vector(3, 0.0), + py::arg("rotation") = Eigen::Matrix3d::Identity(), py::arg("roundrobin") = true, py::arg("posvelrows") = true) .def("initFromArray", [](BasisClasses::BiorthBasis& A, std::vector center) @@ -2112,11 +2117,13 @@ void BasisFactoryClasses(py::module &m) the basis coefficients computed from the particles )", py::arg("reader"), - py::arg("center") = std::vector(3, 0.0)) + py::arg("center") = std::vector(3, 0.0), + py::arg("rotation") = Eigen::Matrix3d::Identity()) .def("initFromArray", - [](BasisClasses::FieldBasis& A, std::vector center) + [](BasisClasses::FieldBasis& A, std::vector center, + Eigen::MatrixXd rotation) { - return A.initFromArray(center); + return A.initFromArray(center, rotation); }, R"( Initialize coefficient accumulation @@ -2125,6 +2132,8 @@ void BasisFactoryClasses(py::module &m) ---------- center : list, default=[0, 0, 0] vector of center positions + rotation : numpy.ndarray, default=Identity + rotation matrix to apply to the phase-space coordinates Returns ------- @@ -2141,7 +2150,8 @@ void BasisFactoryClasses(py::module &m) phase space arrays but allows for generation from very large phase-space sets that can not be stored in physical memory. )", - py::arg("center") = std::vector(3, 0.0)) + py::arg("center") = std::vector(3, 0.0), + py::arg("rotation") = Eigen::Matrix3d::Identity()) .def("addFromArray", [](BasisClasses::FieldBasis& A, Eigen::VectorXd& mass, RowMatrixXd& ps) { diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index c4559acc5..2b1e9edae 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -722,6 +722,11 @@ void CoefficientClasses(py::module &m) { float data's center value )") + .def_readonly("orient", &CoefStruct::rot, + R"( + numpy.ndarray + data's rotation value + )") .def("getCoefTime", &CoefStruct::getTime, R"( Read-only access to the coefficient time @@ -788,6 +793,41 @@ void CoefficientClasses(py::module &m) { Notes ----- + See also + -------- + getCenter : read-only access to center data + )") + .def("getCoefRotation", &CoefStruct::getRotation, + R"( + Read-only access to the rotation matrix + + Returns + ------- + numpy.ndarray + unitary matrix of rotation data + + See also + -------- + setCoefRotation : read-write access to the rotation matrix + )") + .def("setCoefRotation", + static_cast(&CoefStruct::setRotation), + py::arg("mat"), + R"( + Set the rotation matrix + + Parameters + ---------- + mat : numpy.ndarray + rotation matrix + + Returns + ------- + None + + Notes + ----- + See also -------- getCenter : read-only access to center data diff --git a/src/Component.H b/src/Component.H index 8be95eaa0..3dc3b2e8d 100644 --- a/src/Component.H +++ b/src/Component.H @@ -311,7 +311,7 @@ protected: //! Use center from this component Component *c0; - //! Locagte centering component + //! Locate centering component void find_ctr_component(); //! Name of the centering component @@ -766,6 +766,14 @@ public: return ret; } + //! Get rotation matrix + Eigen::Matrix3d getRotation() + { + Eigen::Matrix3d ret = Eigen::Matrix3d::Identity(); + if ( (EJ & Orient::AXIS) and !EJdryrun) ret = orient->transformBody(); + return ret; + } + //! Access to velocities inline double Vel(int i, int j, unsigned flags=Inertial) { PartMap::iterator tp = particles.find(i); From 684f68681e99f19ed00b010ab8442997999a0b04 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Tue, 29 Jul 2025 13:55:55 -0400 Subject: [PATCH 03/18] Update utils/PhaseSpace/KDE2d.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/KDE2d.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/PhaseSpace/KDE2d.cc b/utils/PhaseSpace/KDE2d.cc index 1218d4b1c..9f740858d 100644 --- a/utils/PhaseSpace/KDE2d.cc +++ b/utils/PhaseSpace/KDE2d.cc @@ -114,7 +114,7 @@ namespace KDE for (int j=0; j Date: Tue, 29 Jul 2025 15:14:47 -0400 Subject: [PATCH 04/18] Updated from std::vector to Eigen::Vector3d for center parameter --- expui/BasisFactory.H | 14 ++-- expui/BasisFactory.cc | 2 +- expui/BiorthBasis.H | 6 +- expui/BiorthBasis.cc | 18 ++--- expui/CoefStruct.H | 15 ++-- expui/Coefficients.cc | 46 ++++++----- expui/FieldBasis.H | 4 +- expui/FieldBasis.cc | 4 +- pyEXP/BasisWrappers.cc | 151 ++++++++++++++++++++++++++++++++---- pyEXP/CoefWrappers.cc | 2 +- src/Component.H | 4 +- src/ComponentContainer.cc | 4 +- utils/PhaseSpace/diffpsp.cc | 3 +- 13 files changed, 192 insertions(+), 81 deletions(-) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index c16fffab5..3bb322103 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -66,10 +66,10 @@ namespace BasisClasses unsigned long coefindx; //! The expansion center - std::vector coefctr; + Eigen::Vector3d coefctr = Eigen::Vector3d::Zero(); //! Rotation matrix - Eigen::Matrix3d coefrot; + Eigen::Matrix3d coefrot = Eigen::Matrix3d::Identity(); //! Contains contructor and BFE parameter database YAML::Node node, conf; @@ -154,7 +154,7 @@ namespace BasisClasses virtual ~Basis(void) {} //! Set the expansion center - void setCenter(std::vector center) { coefctr = center; } + void setCenter(Eigen::Vector3d center) { coefctr = center; } //! Set the rotation matrix void setRotation(Eigen::Matrix3d rot) { coefrot = rot; } @@ -218,12 +218,12 @@ namespace BasisClasses //! Generate coeffients from a particle reader virtual CoefClasses::CoefStrPtr - createFromReader(PR::PRptr reader, std::vector ctr, + createFromReader(PR::PRptr reader, Eigen::Vector3d ctr, Eigen::Matrix3d rot) = 0; //! Generate coefficients from a phase-space table virtual void - initFromArray(std::vector ctr, Eigen::Matrix3d rot) = 0; + initFromArray(Eigen::Vector3d ctr, Eigen::Matrix3d rot) = 0; //! Accumulate coefficient contributions from arrays virtual void @@ -234,7 +234,7 @@ namespace BasisClasses //! for the expansion CoefClasses::CoefStrPtr createFromArray (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0, - std::vector center={0.0, 0.0, 0.0}, + Eigen::Vector3d center=Eigen::Vector3d::Zero(), Eigen::Matrix3d rot=Eigen::Matrix3d::Identity(), bool roundrobin=true, bool posvelrows=false); @@ -282,7 +282,7 @@ namespace BasisClasses { return getFieldLabels(coordinates); } //! Get the basis expansion center - std::vector getCenter() { return coefctr; } + Eigen::Vector3d getCenter() { return coefctr; } //! Get the basis expansion center Eigen::Matrix3d getRotation() { return coefrot; } diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index 0e486933b..738ce3f7a 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -275,7 +275,7 @@ namespace BasisClasses // Generate coefficients from a phase-space table // CoefClasses::CoefStrPtr Basis::createFromArray - (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector ctr, + (Eigen::VectorXd& m, RowMatrixXd& p, double time, Eigen::Vector3d ctr, Eigen::Matrix3d rot, bool roundrobin, bool posvelrows) { initFromArray(ctr, rot); diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 7e2f70d12..1c283c8f1 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -94,14 +94,14 @@ namespace BasisClasses //! Generate coeffients from a particle reader and optional center //! location for the expansion CoefClasses::CoefStrPtr createFromReader - (PR::PRptr reader, std::vector center={0.0, 0.0, 0.0}, + (PR::PRptr reader, Eigen::Vector3d center=Eigen::Vector3d::Zero(), Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); //! Generate coeffients from an array and optional center location //! for the expansion CoefClasses::CoefStrPtr createFromArray (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0, - std::vector center={0.0, 0.0, 0.0}, + Eigen::Vector3d center=Eigen::Vector3d::Zero(), Eigen::Matrix3d rot=Eigen::Matrix3d::Identity(), bool roundrobin=true, bool posvelrows=false); @@ -112,7 +112,7 @@ namespace BasisClasses //! Initialize accumulating coefficients from arrays with an optional //! center vector. This is called once to initialize the accumulation. void initFromArray - (std::vector center={0.0, 0.0, 0.0}, + (Eigen::Vector3d center=Eigen::Vector3d::Zero(), Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); //! Initialize accumulating coefficients from arrays. This is diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 5f7a061a4..4d0e78b17 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3871,7 +3871,7 @@ namespace BasisClasses // Generate coeffients from a particle reader CoefClasses::CoefStrPtr BiorthBasis::createFromReader - (PR::PRptr reader, std::vector ctr, Eigen::Matrix3d rot) + (PR::PRptr reader, Eigen::Vector3d ctr, Eigen::Matrix3d rot) { CoefClasses::CoefStrPtr coef; @@ -3922,7 +3922,7 @@ namespace BasisClasses if (use) { Eigen::Vector3d pp; - for (int k=0; k<3; k++) pp(k) = p->pos[k] - coefctr[k]; + for (int k=0; k<3; k++) pp(k) = p->pos[k] - coefctr(k); pp = coefrot * pp; accumulate(pp(0), pp(1), pp(2), p->mass); @@ -3934,7 +3934,7 @@ namespace BasisClasses } // Generate coefficients from a phase-space table - void BiorthBasis::initFromArray(std::vector ctr, Eigen::Matrix3d rot) + void BiorthBasis::initFromArray(Eigen::Vector3d ctr, Eigen::Matrix3d rot) { if (name.compare("sphereSL") == 0) coefret = std::make_shared(); @@ -4038,10 +4038,7 @@ namespace BasisClasses coefindx++; if (use) { - Eigen::Vector3d pp; - for (int k=0; k<3; k++) pp(k) = p(k, n) - coefctr[k]; - pp = coefrot * pp; - + auto pp = coefrot * (p.col(n) - coefctr); accumulate(pp(0), pp(1), pp(2), m(n)); } } @@ -4070,10 +4067,7 @@ namespace BasisClasses coefindx++; if (use) { - Eigen::Vector3d pp; - for (int k=0; k<3; k++) pp(k) = p(k, n) - coefctr[0]; - pp = coefrot * pp; - + auto pp = coefrot * (p.col(n) - coefctr); accumulate(pp(0), pp(1), pp(2), m(n)); } } @@ -4092,7 +4086,7 @@ namespace BasisClasses // Generate coefficients from a phase-space table // CoefClasses::CoefStrPtr BiorthBasis::createFromArray - (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector ctr, + (Eigen::VectorXd& m, RowMatrixXd& p, double time, Eigen::Vector3d ctr, Eigen::Matrix3d rot, bool RoundRobin, bool PosVelRows) { initFromArray(ctr, rot); diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index 052a9e4e0..88e1f2e4b 100644 --- a/expui/CoefStruct.H +++ b/expui/CoefStruct.H @@ -27,7 +27,7 @@ namespace CoefClasses Eigen::VectorXcd store; //! Center data - std::vector ctr= {0.0, 0.0, 0.0}; + Eigen::Vector3d ctr= Eigen::Vector3d::Zero(); //! Orientation data Eigen::Matrix3d rot= Eigen::Matrix3d::Identity(); @@ -65,16 +65,11 @@ namespace CoefClasses //! Read-only access to coefficient data (no copy) Eigen::Ref getCoefs() { return store; } - //! Set new center data (no copy) - void setCenter(std::vector& STORE) - { - if (STORE.size() != ctr.size()) - throw std::invalid_argument("CoefStruct::setCenter: center vector size does not match"); - ctr = STORE; - } + //! Set new center data + void setCenter(Eigen::Vector3d& CTR) {ctr = CTR; } - //! Read-only access to center (no copy) - std::vector getCenter() { return ctr; } + //! Read-only access to center + Eigen::Vector3d getCenter() { return ctr; } //! Set new rotation matrix void setRotation(Eigen::Matrix3d& ROT) { rot = ROT; } diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index 10e2fe233..83403914c 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -113,7 +113,7 @@ namespace CoefClasses // Check for center data // - std::vector ctr; + Eigen::VectorXd ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { stanza.getAttribute("Center").read(ctr); } @@ -148,8 +148,7 @@ namespace CoefClasses // auto coef = std::make_shared(); - if (ctr.size()) coef->ctr = ctr; - + coef->ctr = ctr; coef->rot = rot; coef->lmax = Lmax; coef->nmax = Nmax; @@ -321,7 +320,7 @@ namespace CoefClasses // Check for center data // - std::vector ctr; + Eigen::Vector3d ctr; if (stanza.hasAttribute("Center")) { stanza.getAttribute("Center").read(ctr); } @@ -344,8 +343,7 @@ namespace CoefClasses // auto coef = std::make_shared(); - if (ctr.size()) coef->ctr = ctr; - + coef->ctr = ctr; coef->rot = rot; coef->nfld = Nfld; coef->lmax = Lmax; @@ -420,7 +418,7 @@ namespace CoefClasses // Check for center data // - std::vector ctr; + Eigen::Vector3d ctr; if (stanza.hasAttribute("Center")) { stanza.getAttribute("Center").read(ctr); } @@ -444,8 +442,7 @@ namespace CoefClasses // auto coef = std::make_shared(); - if (ctr.size()) coef->ctr = ctr; - + coef->ctr = ctr; coef->rot = rot; coef->nfld = Nfld; coef->mmax = Mmax; @@ -687,8 +684,7 @@ namespace CoefClasses // Add a center attribute // - if (C->ctr.size()>0) - stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); + stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); // Add a rotation matrix attribute // @@ -870,11 +866,18 @@ namespace CoefClasses // Check for center data // - std::vector ctr; + Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { stanza.getAttribute("Center").read(ctr); } + // Check for rotation data + // + Eigen::Matrix3d rot = Eigen::Matrix3d::Identity(); + if (stanza.hasAttribute("Rotation")) { + stanza.getAttribute("Rotation").read(rot); + } + if (Time < Tmin or Time > Tmax) continue; auto in = stanza.getDataSet("coefficients").read(); @@ -902,8 +905,8 @@ namespace CoefClasses // auto coef = std::make_shared(); - if (ctr.size()) coef->ctr = ctr; - + coef->ctr = ctr; + coef->rot = rot; coef->assign(in, Mmax, Nmax); coef->time = Time; @@ -1083,8 +1086,7 @@ namespace CoefClasses // Add center attribute // - if (C->ctr.size()>0) - stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); + stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); // Add a rotation matrix attribute // @@ -1264,7 +1266,7 @@ namespace CoefClasses // Check for center data // - std::vector ctr; + Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { stanza.getAttribute("Center").read(ctr); } @@ -1616,7 +1618,7 @@ namespace CoefClasses // Check for center data // - std::vector ctr; + Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { stanza.getAttribute("Center").read(ctr); } @@ -2952,10 +2954,8 @@ namespace CoefClasses // Add a center attribute // - if (C->ctr.size()>0) - stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); + stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); - // Add a rotation matrix attribute // Eigen::Matrix3d rot = C->getRotation(); @@ -3077,10 +3077,8 @@ namespace CoefClasses // Add a center attribute // - if (C->ctr.size()>0) - stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); + stanza.createAttribute("Center", HighFive::DataSpace::From(C->ctr)).write(C->ctr); - // Coefficient size (allow Eigen::Tensor to be easily recontructed from metadata) // const auto& d = C->coefs->dimensions(); diff --git a/expui/FieldBasis.H b/expui/FieldBasis.H index 432706c19..efa9877d3 100644 --- a/expui/FieldBasis.H +++ b/expui/FieldBasis.H @@ -135,7 +135,7 @@ namespace BasisClasses //! Generate coeffients from a particle reader and optional center //! location for the expansion CoefClasses::CoefStrPtr createFromReader - (PR::PRptr reader, std::vector center={0.0, 0.0, 0.0}, + (PR::PRptr reader, Eigen::Vector3d center=Eigen::Vector3d::Zero(), Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); @@ -146,7 +146,7 @@ namespace BasisClasses //! Initialize accumulating coefficients from arrays with an optional //! center vector. This is called once to initialize the accumulation. void initFromArray - (std::vector center={0.0, 0.0, 0.0}, + (Eigen::Vector3d center=Eigen::Vector3d::Zero(), Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index a720afaac..6a8f9bf87 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -539,7 +539,7 @@ namespace BasisClasses // Generate coeffients from a particle reader CoefClasses::CoefStrPtr FieldBasis::createFromReader - (PR::PRptr reader, std::vector ctr, Eigen::Matrix3d rot) + (PR::PRptr reader, Eigen::Vector3d ctr, Eigen::Matrix3d rot) { CoefClasses::CoefStrPtr coef; @@ -595,7 +595,7 @@ namespace BasisClasses } // Generate coefficients from a phase-space table - void FieldBasis::initFromArray(std::vector ctr, Eigen::Matrix3d rot) + void FieldBasis::initFromArray(Eigen::Vector3d ctr, Eigen::Matrix3d rot) { if (dof==3) coefret = std::make_shared(); diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 0f91d0f06..58ee68a00 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -292,14 +292,14 @@ void BasisFactoryClasses(py::module &m) } virtual CoefClasses::CoefStrPtr - createFromReader(PR::PRptr reader, std::vector ctr, + createFromReader(PR::PRptr reader, Eigen::Vector3d ctr, Eigen::Matrix3d rot) override { PYBIND11_OVERRIDE_PURE(CoefClasses::CoefStrPtr, Basis, createFromReader, reader, ctr, rot); } virtual void initFromArray - (std::vector ctr, Eigen::MatrixXd rot) { + (Eigen::Vector3d ctr, Eigen::Matrix3d rot) { PYBIND11_OVERRIDE_PURE(void, Basis, initFromArray, ctr, rot); } @@ -888,7 +888,7 @@ void BasisFactoryClasses(py::module &m) ) .def("createFromArray", [](BasisClasses::Basis& A, Eigen::VectorXd& mass, RowMatrixXd& ps, - double time, std::vector center, Eigen::Matrix3d rot, + double time, Eigen::Vector3d center, Eigen::Matrix3d rot, bool roundrobin, bool posvelrows) { return A.createFromArray(mass, ps, time, center, rot, @@ -905,6 +905,11 @@ void BasisFactoryClasses(py::module &m) ps : numpy.ndarray an array with n rows and 6 columns (x, y, z, u, v, w) or 3 columns (x, y, z) for a biorthogonal basis + center : numpy.ndarray + the center of expansion for the basis functions + rotation : numpy.ndarray + the rotation matrix for the basis functions + (default is identity matrix) roundrobin : bool the particles will be accumulated for each process round-robin style with MPI by default. This may be @@ -927,9 +932,53 @@ void BasisFactoryClasses(py::module &m) makeFromArray : create coefficients contributions )", py::arg("mass"), py::arg("pos"), py::arg("time"), - py::arg("center") = std::vector(3, 0.0), + py::arg("center") = Eigen::Vector3d::Zero(), py::arg("rotation") = Eigen::Matrix3d::Identity(), py::arg("roundrobin") = true, py::arg("posvelrows") = true) + .def("createFromArray", + [](BasisClasses::Basis& A, Eigen::VectorXd& mass, RowMatrixXd& ps, + double time, Eigen::Vector3d center, bool roundrobin, bool posvelrows) + { + return A.createFromArray(mass, ps, time, center, Eigen::Matrix3d::Identity(), + roundrobin, posvelrows); + }, + R"( + Generate the coefficients from a mass and position array or, + phase-space array, time, and an optional expansion center location. + + Parameters + ---------- + mass : list + vector containing the masses for the n particles + ps : numpy.ndarray + an array with n rows and 6 columns (x, y, z, u, v, w) + or 3 columns (x, y, z) for a biorthogonal basis + center : numpy.ndarray + the center of expansion for the basis functions + roundrobin : bool + the particles will be accumulated for each process + round-robin style with MPI by default. This may be + disabled with 'roundrobin=False'. + posvelrows : bool + positions (and optionally velocities) will be packed + in rows instead of columns. This accommodates the numpy + construction [xpos, ypos, zpos] where xpos, ypos, zpos are + arrays. Defaults to True. + + Returns + ------- + CoefStruct + the coefficient structure derived from the particles + + See also + -------- + initFromArray : initialize for coefficient contributions + addFromArray : add contribution for particles + makeFromArray : create coefficients contributions + )", + py::arg("mass"), py::arg("pos"), py::arg("time"), + py::arg("center") = Eigen::Vector3d::Zero(), + py::arg("roundrobin") = true, py::arg("posvelrows") = true) .def("makeFromArray", [](BasisClasses::Basis& A, double time) { @@ -1130,11 +1179,11 @@ void BasisFactoryClasses(py::module &m) the basis coefficients computed from the particles )", py::arg("reader"), - py::arg("center") = std::vector(3, 0.0), + py::arg("center") = Eigen::Vector3d::Zero(), py::arg("rotation") = Eigen::Matrix3d::Identity()) .def("createFromArray", [](BasisClasses::BiorthBasis& A, Eigen::VectorXd& mass, RowMatrixXd& pos, - double time, std::vector center, Eigen::Matrix3d rot, + double time, Eigen::Vector3d center, Eigen::Matrix3d rot, bool roundrobin, bool posvelrows) { return A.createFromArray(mass, pos, time, center, rot, @@ -1172,13 +1221,57 @@ void BasisFactoryClasses(py::module &m) makeFromArray : create coefficients contributions )", py::arg("mass"), py::arg("pos"), py::arg("time"), - py::arg("center") = std::vector(3, 0.0), + py::arg("center") = Eigen::Vector3d::Zero(), py::arg("rotation") = Eigen::Matrix3d::Identity(), py::arg("roundrobin") = true, py::arg("posvelrows") = true) + .def("createFromArray", + [](BasisClasses::BiorthBasis& A, Eigen::VectorXd& mass, RowMatrixXd& pos, + double time, Eigen::Vector3d center, bool roundrobin, bool posvelrows) + { + return A.createFromArray(mass, pos, time, center, Eigen::Matrix3d::Identity(), + roundrobin, posvelrows); + }, + R"( + Generate the coefficients from a mass and position array, + time, and an optional expansion center location. + + Parameters + ---------- + mass : list + vector containing the masses for the n particles + pos : numpy.ndarray + an array with n rows and 3 columns (x, y, z) + center : numpy.ndarray + the center of expansion for the basis functions + roundrobin : bool + the particles will be accumulated for each process + round-robin style with MPI by default. This may be + disabled with 'roundrobin=false'. + posvelrows : bool + positions (and optionally velocities) will be packed + in rows instead of columns. This accommodates the numpy + construction [xpos, ypos, zpos] where xpos, ypos, zpos are + arrays. Defaults to True. + + Returns + ------- + CoefStruct + the coefficient structure derived from the particles + + See also + -------- + initFromArray : initialize for coefficient contributions + addFromArray : add contribution for particles + makeFromArray : create coefficients contributions + )", + py::arg("mass"), py::arg("pos"), py::arg("time"), + py::arg("center") = Eigen::Vector3d::Zero(), + py::arg("roundrobin") = true, py::arg("posvelrows") = true) .def("initFromArray", - [](BasisClasses::BiorthBasis& A, std::vector center) + [](BasisClasses::BiorthBasis& A, + Eigen::Vector3d ctr, Eigen::Matrix3d rot) { - return A.initFromArray(center); + return A.initFromArray(ctr, rot); }, R"( Initialize coefficient accumulation @@ -1203,7 +1296,8 @@ void BasisFactoryClasses(py::module &m) phase space arrays but allows for generation from very large phase-space sets that can not be stored in physical memory. )", - py::arg("center") = std::vector(3, 0.0)) + py::arg("center") = Eigen::Vector3d::Zero(), + py::arg("rotation") = Eigen::Matrix3d::Identity()) .def("addFromArray", [](BasisClasses::BiorthBasis& A, Eigen::VectorXd& mass, RowMatrixXd& pos) { @@ -2120,8 +2214,8 @@ void BasisFactoryClasses(py::module &m) py::arg("center") = std::vector(3, 0.0), py::arg("rotation") = Eigen::Matrix3d::Identity()) .def("initFromArray", - [](BasisClasses::FieldBasis& A, std::vector center, - Eigen::MatrixXd rotation) + [](BasisClasses::FieldBasis& A, + Eigen::Vector3d center, Eigen::MatrixXd rotation) { return A.initFromArray(center, rotation); }, @@ -2150,8 +2244,39 @@ void BasisFactoryClasses(py::module &m) phase space arrays but allows for generation from very large phase-space sets that can not be stored in physical memory. )", - py::arg("center") = std::vector(3, 0.0), + py::arg("center") = Eigen::Vector3d::Zero(), py::arg("rotation") = Eigen::Matrix3d::Identity()) + .def("initFromArray", + [](BasisClasses::FieldBasis& A, Eigen::Vector3d center) + { + return A.initFromArray(center, Eigen::Matrix3d::Identity()); + }, + R"( + Initialize coefficient accumulation + + Parameters + ---------- + center : list, default=[0, 0, 0] + vector of center positions + rotation : numpy.ndarray, default=Identity + rotation matrix to apply to the phase-space coordinates + + Returns + ------- + None + + Notes + ----- + After initialization, phase-space data is then added with + addFromArray() call. addFromArray() may be called multiple times + with any unique partition of phase space. The final generation is + finished with a call to makeFromArray() with the snapshot time. + This final call returns the coefficient set. This sequence of + calls is identical to createFromArray() for a single set of + phase space arrays but allows for generation from very large + phase-space sets that can not be stored in physical memory. + )", + py::arg("center") = Eigen::Vector3d::Zero()) .def("addFromArray", [](BasisClasses::FieldBasis& A, Eigen::VectorXd& mass, RowMatrixXd& ps) { diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index 2b1e9edae..48d344049 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -776,7 +776,7 @@ void CoefficientClasses(py::module &m) { setCenter : read-write access to the center data )") .def("setCoefCenter", - static_cast&)>(&CoefStruct::setCenter), + static_cast(&CoefStruct::setCenter), py::arg("mat"), R"( Set the center vector diff --git a/src/Component.H b/src/Component.H index 3dc3b2e8d..0edcea661 100644 --- a/src/Component.H +++ b/src/Component.H @@ -755,9 +755,9 @@ public: } //! Get rectified component center - std::vector getCenter(unsigned flags=Inertial) + Eigen::Vector3d getCenter(unsigned flags=Inertial) { - std::vector ret(3, 0.0); + Eigen::Vector3d ret = Eigen::Vector3d::Zero(); for (int j=0; j<3; j++) { if (com_system and flags & Local) ret[j] += com0[j]; if (flags & Centered) ret[j] += center[j]; diff --git a/src/ComponentContainer.cc b/src/ComponentContainer.cc index 8f95a49d3..d64885e0c 100644 --- a/src/ComponentContainer.cc +++ b/src/ComponentContainer.cc @@ -823,8 +823,8 @@ void ComponentContainer::compute_potential(unsigned mlevel) other->time_so_far.stop(); if (false) { // Some deep debugging for playback . . . - std::vector cen1 = inter->c->getCenter(Component::Local); - std::vector cen2 = other->getCenter(Component::Local); + auto cen1 = inter->c->getCenter(Component::Local); + auto cen2 = other->getCenter(Component::Local); std::cout << "ComponentContainer [" << myid << "], centers for [" << inter->c->name << "-->" << other->name << "] c1=(" diff --git a/utils/PhaseSpace/diffpsp.cc b/utils/PhaseSpace/diffpsp.cc index 7a52b03bd..c7b2a5427 100755 --- a/utils/PhaseSpace/diffpsp.cc +++ b/utils/PhaseSpace/diffpsp.cc @@ -630,10 +630,9 @@ main(int argc, char **argv) // Number of paths // - int npath1 = 1, npath2 = 1; + int npath1 = 1; if (fileType != "PSPhdf5") { npath1 = INFILE1.size(); - npath2 = INFILE1.size(); } // Iterate through file list From 3f6e5aa48b4323b411a1c3ccbc7936095980af47 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 29 Jul 2025 16:18:30 -0400 Subject: [PATCH 05/18] Exploit scalar*vector algebra in Eigen --- expui/BiorthBasis.cc | 22 ++++++---------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 4d0e78b17..e82a994db 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4205,28 +4205,18 @@ namespace BasisClasses // Interpolate center // - if (coefsA->ctr.size() and coefsB->ctr.size()) { - newcoef->ctr.resize(3); - for (int k=0; k<3; k++) - newcoef->ctr[k] = a * coefsA->ctr[k] + b * coefsB->ctr[k]; - } + newcoef->ctr = a * coefsA->ctr + b * coefsB->ctr; // Interpolate rotation matrix followed by unitarization // - Eigen::Matrix3d newrot; - for (int k=0; k<9; k++) { - newrot.data()[k] = - a * coefsA->rot.data()[k] + b * coefsB->rot.data()[k]; - } - + Eigen::Matrix3d newrot = a * coefsA->rot + b * coefsB->rot; + // Closest unitary matrix in the Frobenius norm sense // - Eigen::BDCSVD svd(newrot, - Eigen::ComputeFullU | Eigen::ComputeFullV); - auto U = svd.matrixU(); - auto V = svd.matrixV(); + Eigen::BDCSVD svd + (newrot, Eigen::ComputeFullU | Eigen::ComputeFullV); - newcoef->rot = U * V.adjoint(); + newcoef->rot = svd.matrixU() * svd.matrixV().adjoint(); // Install coefficients // From a2547bcca6b85f8dc4764871e2af94fc0fdc3659 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 30 Jul 2025 10:23:08 -0400 Subject: [PATCH 06/18] Work around for HighFive std::vector and Eigen::Vector dimensionality in attribute reading/writing --- expui/Coefficients.cc | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index 83403914c..1fba2b738 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -115,7 +115,9 @@ namespace CoefClasses // Eigen::VectorXd ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { - stanza.getAttribute("Center").read(ctr); + std::vector ctrvec; + stanza.getAttribute("Center").read(ctrvec); + for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; } // Check for rotation matrix @@ -322,7 +324,9 @@ namespace CoefClasses // Eigen::Vector3d ctr; if (stanza.hasAttribute("Center")) { + std::vector ctrvec; stanza.getAttribute("Center").read(ctr); + for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; } // Check for rotation matrix @@ -420,7 +424,9 @@ namespace CoefClasses // Eigen::Vector3d ctr; if (stanza.hasAttribute("Center")) { + std::vector ctrvec; stanza.getAttribute("Center").read(ctr); + for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; } // Check for rotation matrix @@ -868,7 +874,9 @@ namespace CoefClasses // Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { + std::vector ctrvec; stanza.getAttribute("Center").read(ctr); + for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; } // Check for rotation data @@ -1268,7 +1276,9 @@ namespace CoefClasses // Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { + std::vector ctrvec; stanza.getAttribute("Center").read(ctr); + for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; } if (Time < Tmin or Time > Tmax) continue; @@ -1620,7 +1630,9 @@ namespace CoefClasses // Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { + std::vector ctrvec; stanza.getAttribute("Center").read(ctr); + for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; } if (Time < Tmin or Time > Tmax) continue; From 15ecc684269f96ed48593bbacd90eb30f3cee720 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 30 Jul 2025 13:12:03 -0400 Subject: [PATCH 07/18] Added a comment on the Center autodeduction workaround for HighFive --- expui/Coefficients.cc | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index 1fba2b738..78ff0b3d9 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -115,6 +115,10 @@ namespace CoefClasses // Eigen::VectorXd ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { + // The previous version used a std::vector for the + // center which HighFive autodeduced to a 1D data space, so we + // need to read it with a std::vector for backwards + // compatibility std::vector ctrvec; stanza.getAttribute("Center").read(ctrvec); for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; @@ -324,6 +328,10 @@ namespace CoefClasses // Eigen::Vector3d ctr; if (stanza.hasAttribute("Center")) { + // The previous version used a std::vector for the + // center which HighFive autodeduced to a 1D data space, so we + // need to read it with a std::vector for backwards + // compatibility std::vector ctrvec; stanza.getAttribute("Center").read(ctr); for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; @@ -424,6 +432,10 @@ namespace CoefClasses // Eigen::Vector3d ctr; if (stanza.hasAttribute("Center")) { + // The previous version used a std::vector for the + // center which HighFive autodeduced to a 1D data space, so we + // need to read it with a std::vector for backwards + // compatibility std::vector ctrvec; stanza.getAttribute("Center").read(ctr); for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; @@ -874,6 +886,10 @@ namespace CoefClasses // Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { + // The previous version used a std::vector for the + // center which HighFive autodeduced to a 1D data space, so we + // need to read it with a std::vector for backwards + // compatibility std::vector ctrvec; stanza.getAttribute("Center").read(ctr); for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; @@ -1276,6 +1292,10 @@ namespace CoefClasses // Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { + // The previous version used a std::vector for the + // center which HighFive autodeduced to a 1D data space, so we + // need to read it with a std::vector for backwards + // compatibility std::vector ctrvec; stanza.getAttribute("Center").read(ctr); for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; @@ -1630,6 +1650,9 @@ namespace CoefClasses // Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { + // The previous version used a std::vector for the + // center which autodeduced to a 1D vector, so we need to read + // it with a std::vector for backwards compatibility std::vector ctrvec; stanza.getAttribute("Center").read(ctr); for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; From 6d876e765253d5326aa0143ce9d20eec87981a57 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 31 Jul 2025 11:13:57 -0400 Subject: [PATCH 08/18] Ignore dimensionality when reading center vector --- expui/Coefficients.cc | 59 +++++++++++-------------------------------- 1 file changed, 15 insertions(+), 44 deletions(-) diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index 78ff0b3d9..86fd19e5a 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -113,15 +113,10 @@ namespace CoefClasses // Check for center data // - Eigen::VectorXd ctr = Eigen::Vector3d::Zero(); + Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { - // The previous version used a std::vector for the - // center which HighFive autodeduced to a 1D data space, so we - // need to read it with a std::vector for backwards - // compatibility - std::vector ctrvec; - stanza.getAttribute("Center").read(ctrvec); - for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; + // Read three values, ignoring dimensionality + stanza.getAttribute("Center").read_raw(ctr.data()); } // Check for rotation matrix @@ -326,15 +321,10 @@ namespace CoefClasses // Check for center data // - Eigen::Vector3d ctr; + Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { - // The previous version used a std::vector for the - // center which HighFive autodeduced to a 1D data space, so we - // need to read it with a std::vector for backwards - // compatibility - std::vector ctrvec; - stanza.getAttribute("Center").read(ctr); - for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; + // Read three values, ignoring dimensionality + stanza.getAttribute("Center").read_raw(ctr.data()); } // Check for rotation matrix @@ -430,15 +420,10 @@ namespace CoefClasses // Check for center data // - Eigen::Vector3d ctr; + Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { - // The previous version used a std::vector for the - // center which HighFive autodeduced to a 1D data space, so we - // need to read it with a std::vector for backwards - // compatibility - std::vector ctrvec; - stanza.getAttribute("Center").read(ctr); - for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; + // Read three values, ignoring dimensionality + stanza.getAttribute("Center").read_raw(ctr.data()); } // Check for rotation matrix @@ -886,13 +871,8 @@ namespace CoefClasses // Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { - // The previous version used a std::vector for the - // center which HighFive autodeduced to a 1D data space, so we - // need to read it with a std::vector for backwards - // compatibility - std::vector ctrvec; - stanza.getAttribute("Center").read(ctr); - for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; + // Read three values, ignoring dimensionality + stanza.getAttribute("Center").read_raw(ctr.data()); } // Check for rotation data @@ -1292,13 +1272,8 @@ namespace CoefClasses // Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { - // The previous version used a std::vector for the - // center which HighFive autodeduced to a 1D data space, so we - // need to read it with a std::vector for backwards - // compatibility - std::vector ctrvec; - stanza.getAttribute("Center").read(ctr); - for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; + // Read three values, ignoring dimensionality + stanza.getAttribute("Center").read_raw(ctr.data()); } if (Time < Tmin or Time > Tmax) continue; @@ -1650,12 +1625,8 @@ namespace CoefClasses // Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); if (stanza.hasAttribute("Center")) { - // The previous version used a std::vector for the - // center which autodeduced to a 1D vector, so we need to read - // it with a std::vector for backwards compatibility - std::vector ctrvec; - stanza.getAttribute("Center").read(ctr); - for (int k=0; k<3; k++) ctr(k) = ctrvec[k]; + // Read three values, ignoring dimensionality + stanza.getAttribute("Center").read_raw(ctr.data()); } if (Time < Tmin or Time > Tmax) continue; From 6913257e9eb685c0aed7df42d82cbbfd917aa1b0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 6 Aug 2025 14:14:38 -0400 Subject: [PATCH 09/18] Wrong type in HighFive attribute creation --- expui/Coefficients.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index 86fd19e5a..0bb08943c 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -1095,7 +1095,7 @@ namespace CoefClasses // Add a rotation matrix attribute // Eigen::Matrix3d rot = C->getRotation(); - stanza.createAttribute("Rotation", HighFive::DataSpace::From(rot)).write(rot); + stanza.createAttribute("Rotation", HighFive::DataSpace::From(rot)).write(rot); // Add coefficient data // @@ -2965,7 +2965,7 @@ namespace CoefClasses // Add a rotation matrix attribute // Eigen::Matrix3d rot = C->getRotation(); - stanza.createAttribute("Rotation", HighFive::DataSpace::From(rot)).write(rot); + stanza.createAttribute("Rotation", HighFive::DataSpace::From(rot)).write(rot); // Coefficient size (allow Eigen::Tensor to be easily recontructed from metadata) // From cdb5d78b23c9683d070ee1a0067f7e1e7f46fc8c Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Sat, 16 Aug 2025 09:49:14 +0100 Subject: [PATCH 10/18] Update pyEXP/BasisWrappers.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 58ee68a00..186e27f02 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2211,7 +2211,7 @@ void BasisFactoryClasses(py::module &m) the basis coefficients computed from the particles )", py::arg("reader"), - py::arg("center") = std::vector(3, 0.0), + py::arg("center") = Eigen::Vector3d::Zero(), py::arg("rotation") = Eigen::Matrix3d::Identity()) .def("initFromArray", [](BasisClasses::FieldBasis& A, From 34a5692a8fd33c0b9d53235caff13dd0975b3e64 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Sat, 16 Aug 2025 09:49:49 +0100 Subject: [PATCH 11/18] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index e82a994db..eb54c7caf 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4116,7 +4116,7 @@ namespace BasisClasses int rows = accel.rows(); for (int n=0; ngetFields(pp(0), pp(1), pp(2)); From a1ce38b0e7d728ed65d300f3cffc213efeefe7ea Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 16 Aug 2025 10:58:02 -0400 Subject: [PATCH 12/18] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index eb54c7caf..8d84eaad8 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4067,7 +4067,7 @@ namespace BasisClasses coefindx++; if (use) { - auto pp = coefrot * (p.col(n) - coefctr); + auto pp = coefrot * (p.row(n).transpose() - coefctr); accumulate(pp(0), pp(1), pp(2), m(n)); } } From 2dadc179293af263acffdf0d879ccf495355f7ee Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 16 Aug 2025 16:15:56 -0400 Subject: [PATCH 13/18] Another fix to the posvel input logic --- expui/BiorthBasis.cc | 60 +++++++++++++++++++++++++++++++++----------- 1 file changed, 46 insertions(+), 14 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 8d84eaad8..761d72f35 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3905,26 +3905,35 @@ namespace BasisClasses // coef->rot = rot; - std::vector pp(3), vv(3); + std::vector p1(3), v1(3); + + // Map the vector rather than copy + // + Eigen::Map pp(p1.data(), 3), vv(v1.data(), 3); + vv.setZero(); reset_coefs(); for (auto p=reader->firstParticle(); p!=0; p=reader->nextParticle()) { bool use = false; + // Translate and rotate the position vector + // + for (int k=0; k<3; k++) p1[k] = p->pos[k]; + pp = coefrot * (pp - coefctr); + if (ftor) { - pp.assign(p->pos, p->pos+3); - vv.assign(p->vel, p->vel+3); - use = ftor(p->mass, pp, vv, p->indx); + // Rotate the velocity vector + // + for (int k=0; k<3; k++) v1[k] = p->vel[k]; + vv = coefrot * vv; + + use = ftor(p->mass, p1, v1, p->indx); } else { use = true; } if (use) { - Eigen::Vector3d pp; - for (int k=0; k<3; k++) pp(k) = p->pos[k] - coefctr(k); - pp = coefrot * pp; - accumulate(pp(0), pp(1), pp(2), p->mass); } } @@ -3994,6 +4003,7 @@ namespace BasisClasses int cols = p.cols(); bool ambiguous = false; + bool haveVel = false; if (cols==3 or cols==6) { if (rows != 3 and rows != 6) PosVelRows = false; @@ -4014,7 +4024,11 @@ namespace BasisClasses << "if this assumption is wrong." << std::endl; } - std::vector p1(3), v1(3, 0); + // Map the vector rather than copy + // + std::vector p1(3), v1(3); + Eigen::Map pp(p1.data(), 3), vv(v1.data(), 3); + vv.setZero(); if (PosVelRows) { if (p.rows()<3) { @@ -4024,13 +4038,23 @@ namespace BasisClasses throw std::runtime_error(msg.str()); } - for (int n=0; n Date: Sun, 17 Aug 2025 19:33:54 -0400 Subject: [PATCH 14/18] Fix a missing registration of the rotation matrix; interface with numpy in row-major space --- expui/BasisFactory.H | 14 ++++--- expui/BasisFactory.cc | 2 +- expui/BiorthBasis.H | 6 +-- expui/BiorthBasis.cc | 45 +++++++++------------- expui/FieldBasis.H | 4 +- expui/FieldBasis.cc | 4 +- pyEXP/BasisWrappers.cc | 84 ++++++++++++++++-------------------------- 7 files changed, 64 insertions(+), 95 deletions(-) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index 3bb322103..b0480da39 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -24,6 +24,8 @@ namespace BasisClasses { using RowMatrixXd = Eigen::Matrix; + using RowMatrix3d = Eigen::Matrix; + //! Callback function signature for selection particles to //! accumulate using Callback = @@ -69,7 +71,7 @@ namespace BasisClasses Eigen::Vector3d coefctr = Eigen::Vector3d::Zero(); //! Rotation matrix - Eigen::Matrix3d coefrot = Eigen::Matrix3d::Identity(); + RowMatrix3d coefrot = Eigen::Matrix3d::Identity(); //! Contains contructor and BFE parameter database YAML::Node node, conf; @@ -157,7 +159,7 @@ namespace BasisClasses void setCenter(Eigen::Vector3d center) { coefctr = center; } //! Set the rotation matrix - void setRotation(Eigen::Matrix3d rot) { coefrot = rot; } + void setRotation(RowMatrix3d rot) { coefrot = rot; } //! Evaluate basis in desired coordinates virtual std::vector @@ -219,11 +221,11 @@ namespace BasisClasses //! Generate coeffients from a particle reader virtual CoefClasses::CoefStrPtr createFromReader(PR::PRptr reader, Eigen::Vector3d ctr, - Eigen::Matrix3d rot) = 0; + RowMatrix3d rot) = 0; //! Generate coefficients from a phase-space table virtual void - initFromArray(Eigen::Vector3d ctr, Eigen::Matrix3d rot) = 0; + initFromArray(Eigen::Vector3d ctr, RowMatrix3d rot) = 0; //! Accumulate coefficient contributions from arrays virtual void @@ -235,7 +237,7 @@ namespace BasisClasses CoefClasses::CoefStrPtr createFromArray (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0, Eigen::Vector3d center=Eigen::Vector3d::Zero(), - Eigen::Matrix3d rot=Eigen::Matrix3d::Identity(), + RowMatrix3d rot=RowMatrix3d::Identity(), bool roundrobin=true, bool posvelrows=false); //! Create and the coefficients from the array accumulation with the @@ -285,7 +287,7 @@ namespace BasisClasses Eigen::Vector3d getCenter() { return coefctr; } //! Get the basis expansion center - Eigen::Matrix3d getRotation() { return coefrot; } + RowMatrix3d getRotation() { return coefrot; } }; using BasisPtr = std::shared_ptr; diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index 738ce3f7a..49b5fe728 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -276,7 +276,7 @@ namespace BasisClasses // CoefClasses::CoefStrPtr Basis::createFromArray (Eigen::VectorXd& m, RowMatrixXd& p, double time, Eigen::Vector3d ctr, - Eigen::Matrix3d rot, bool roundrobin, bool posvelrows) + RowMatrix3d rot, bool roundrobin, bool posvelrows) { initFromArray(ctr, rot); addFromArray(m, p, roundrobin, posvelrows); diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 1c283c8f1..4d40e6b5d 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -95,14 +95,14 @@ namespace BasisClasses //! location for the expansion CoefClasses::CoefStrPtr createFromReader (PR::PRptr reader, Eigen::Vector3d center=Eigen::Vector3d::Zero(), - Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); + RowMatrix3d rot=RowMatrix3d::Identity()); //! Generate coeffients from an array and optional center location //! for the expansion CoefClasses::CoefStrPtr createFromArray (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0, Eigen::Vector3d center=Eigen::Vector3d::Zero(), - Eigen::Matrix3d rot=Eigen::Matrix3d::Identity(), + RowMatrix3d rot=RowMatrix3d::Identity(), bool roundrobin=true, bool posvelrows=false); //! Generate coeffients from an array and optional center location @@ -113,7 +113,7 @@ namespace BasisClasses //! center vector. This is called once to initialize the accumulation. void initFromArray (Eigen::Vector3d center=Eigen::Vector3d::Zero(), - Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); + RowMatrix3d rot=RowMatrix3d::Identity()); //! Initialize accumulating coefficients from arrays. This is //! called once to initialize the accumulation. diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 761d72f35..b2b507717 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3871,7 +3871,7 @@ namespace BasisClasses // Generate coeffients from a particle reader CoefClasses::CoefStrPtr BiorthBasis::createFromReader - (PR::PRptr reader, Eigen::Vector3d ctr, Eigen::Matrix3d rot) + (PR::PRptr reader, Eigen::Vector3d ctr, RowMatrix3d rot) { CoefClasses::CoefStrPtr coef; @@ -3890,19 +3890,14 @@ namespace BasisClasses throw std::runtime_error(sout.str()); } - // Is center non-zero? + // Add the expansion center metadata and register for this instance // - bool addCenter = false; - for (auto v : ctr) { - if (v != 0.0) addCenter = true; - } - - // Add the expansion center metadata - // - if (addCenter) coef->ctr = ctr; + coefctr = ctr; + coef->ctr = ctr; - // Add the rotation matrix + // Add the rotation matrix metadata and register for this instance // + coefrot = rot; coef->rot = rot; std::vector p1(3), v1(3); @@ -3943,7 +3938,7 @@ namespace BasisClasses } // Generate coefficients from a phase-space table - void BiorthBasis::initFromArray(Eigen::Vector3d ctr, Eigen::Matrix3d rot) + void BiorthBasis::initFromArray(Eigen::Vector3d ctr, RowMatrix3d rot) { if (name.compare("sphereSL") == 0) coefret = std::make_shared(); @@ -3958,22 +3953,13 @@ namespace BasisClasses throw std::runtime_error(sout.str()); } - // Is center non-zero? - // - bool addCenter = false; - for (auto v : ctr) { - if (v != 0.0) addCenter = true; - } - - // Add the expansion center metadata - // - if (addCenter) coefret->ctr = ctr; - - // Register the center + // Add the expansion center metadata and register // coefctr = ctr; + coefret->ctr = ctr; - // Register the rotation matrix + // Add the rotation metadata and register + coefrot = rot; coefret->rot = rot; // Clean up for accumulation @@ -4030,6 +4016,9 @@ namespace BasisClasses Eigen::Map pp(p1.data(), 3), vv(v1.data(), 3); vv.setZero(); + if (myid==0) std::cout << "Center: " << coefctr << std::endl + << "Orient: " << coefrot << std::endl; + if (PosVelRows) { if (p.rows()<3) { std::ostringstream msg; @@ -4119,7 +4108,7 @@ namespace BasisClasses // CoefClasses::CoefStrPtr BiorthBasis::createFromArray (Eigen::VectorXd& m, RowMatrixXd& p, double time, Eigen::Vector3d ctr, - Eigen::Matrix3d rot, bool RoundRobin, bool PosVelRows) + RowMatrix3d rot, bool RoundRobin, bool PosVelRows) { initFromArray(ctr, rot); addFromArray(m, p, RoundRobin, PosVelRows); @@ -4241,11 +4230,11 @@ namespace BasisClasses // Interpolate rotation matrix followed by unitarization // - Eigen::Matrix3d newrot = a * coefsA->rot + b * coefsB->rot; + RowMatrix3d newrot = a * coefsA->rot + b * coefsB->rot; // Closest unitary matrix in the Frobenius norm sense // - Eigen::BDCSVD svd + Eigen::BDCSVD svd (newrot, Eigen::ComputeFullU | Eigen::ComputeFullV); newcoef->rot = svd.matrixU() * svd.matrixV().adjoint(); diff --git a/expui/FieldBasis.H b/expui/FieldBasis.H index efa9877d3..37fc20044 100644 --- a/expui/FieldBasis.H +++ b/expui/FieldBasis.H @@ -136,7 +136,7 @@ namespace BasisClasses //! location for the expansion CoefClasses::CoefStrPtr createFromReader (PR::PRptr reader, Eigen::Vector3d center=Eigen::Vector3d::Zero(), - Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); + RowMatrix3d rot=RowMatrix3d::Identity()); //! Generate coeffients from an array and optional center location @@ -147,7 +147,7 @@ namespace BasisClasses //! center vector. This is called once to initialize the accumulation. void initFromArray (Eigen::Vector3d center=Eigen::Vector3d::Zero(), - Eigen::Matrix3d rot=Eigen::Matrix3d::Identity()); + RowMatrix3d rot=RowMatrix3d::Identity()); //! Initialize accumulating coefficients from arrays. This is diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 6a8f9bf87..66a928f34 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -539,7 +539,7 @@ namespace BasisClasses // Generate coeffients from a particle reader CoefClasses::CoefStrPtr FieldBasis::createFromReader - (PR::PRptr reader, Eigen::Vector3d ctr, Eigen::Matrix3d rot) + (PR::PRptr reader, Eigen::Vector3d ctr, RowMatrix3d rot) { CoefClasses::CoefStrPtr coef; @@ -595,7 +595,7 @@ namespace BasisClasses } // Generate coefficients from a phase-space table - void FieldBasis::initFromArray(Eigen::Vector3d ctr, Eigen::Matrix3d rot) + void FieldBasis::initFromArray(Eigen::Vector3d ctr, RowMatrix3d rot) { if (dof==3) coefret = std::make_shared(); diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 186e27f02..f2fd10a5a 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -293,13 +293,13 @@ void BasisFactoryClasses(py::module &m) virtual CoefClasses::CoefStrPtr createFromReader(PR::PRptr reader, Eigen::Vector3d ctr, - Eigen::Matrix3d rot) override { + RowMatrix3d rot) override { PYBIND11_OVERRIDE_PURE(CoefClasses::CoefStrPtr, Basis, createFromReader, reader, ctr, rot); } virtual void initFromArray - (Eigen::Vector3d ctr, Eigen::Matrix3d rot) { + (Eigen::Vector3d ctr, RowMatrix3d rot) { PYBIND11_OVERRIDE_PURE(void, Basis, initFromArray, ctr, rot); } @@ -886,14 +886,7 @@ void BasisFactoryClasses(py::module &m) This is an experimental feature )" ) - .def("createFromArray", - [](BasisClasses::Basis& A, Eigen::VectorXd& mass, RowMatrixXd& ps, - double time, Eigen::Vector3d center, Eigen::Matrix3d rot, - bool roundrobin, bool posvelrows) - { - return A.createFromArray(mass, ps, time, center, rot, - roundrobin, posvelrows); - }, + .def("createFromArray", &BasisClasses::Basis::createFromArray, R"( Generate the coefficients from a mass and position array or, phase-space array, time, and an optional expansion center location. @@ -933,13 +926,13 @@ void BasisFactoryClasses(py::module &m) )", py::arg("mass"), py::arg("pos"), py::arg("time"), py::arg("center") = Eigen::Vector3d::Zero(), - py::arg("rotation") = Eigen::Matrix3d::Identity(), + py::arg("rotation") = RowMatrix3d::Identity(), py::arg("roundrobin") = true, py::arg("posvelrows") = true) .def("createFromArray", [](BasisClasses::Basis& A, Eigen::VectorXd& mass, RowMatrixXd& ps, double time, Eigen::Vector3d center, bool roundrobin, bool posvelrows) { - return A.createFromArray(mass, ps, time, center, Eigen::Matrix3d::Identity(), + return A.createFromArray(mass, ps, time, center, RowMatrix3d::Identity(), roundrobin, posvelrows); }, R"( @@ -979,11 +972,7 @@ void BasisFactoryClasses(py::module &m) py::arg("mass"), py::arg("pos"), py::arg("time"), py::arg("center") = Eigen::Vector3d::Zero(), py::arg("roundrobin") = true, py::arg("posvelrows") = true) - .def("makeFromArray", - [](BasisClasses::Basis& A, double time) - { - return A.makeFromArray(time); - }, + .def("makeFromArray", &BasisClasses::Basis::makeFromArray, R"( Make the coefficients @@ -1172,6 +1161,8 @@ void BasisFactoryClasses(py::module &m) the ParticleReader instance center : list, default=[0, 0, 0] an optional expansion center location + rotation : numpy.ndarray, default=Identity + rotation matrix to apply to the phase-space coordinates Returns ------- @@ -1180,15 +1171,8 @@ void BasisFactoryClasses(py::module &m) )", py::arg("reader"), py::arg("center") = Eigen::Vector3d::Zero(), - py::arg("rotation") = Eigen::Matrix3d::Identity()) - .def("createFromArray", - [](BasisClasses::BiorthBasis& A, Eigen::VectorXd& mass, RowMatrixXd& pos, - double time, Eigen::Vector3d center, Eigen::Matrix3d rot, - bool roundrobin, bool posvelrows) - { - return A.createFromArray(mass, pos, time, center, rot, - roundrobin, posvelrows); - }, + py::arg("rotation") = RowMatrix3d::Identity()) + .def("createFromArray", &BasisClasses::BiorthBasis::createFromArray, R"( Generate the coefficients from a mass and position array, time, and an optional expansion center location. @@ -1222,13 +1206,13 @@ void BasisFactoryClasses(py::module &m) )", py::arg("mass"), py::arg("pos"), py::arg("time"), py::arg("center") = Eigen::Vector3d::Zero(), - py::arg("rotation") = Eigen::Matrix3d::Identity(), + py::arg("rotation") = RowMatrix3d::Identity(), py::arg("roundrobin") = true, py::arg("posvelrows") = true) .def("createFromArray", [](BasisClasses::BiorthBasis& A, Eigen::VectorXd& mass, RowMatrixXd& pos, double time, Eigen::Vector3d center, bool roundrobin, bool posvelrows) { - return A.createFromArray(mass, pos, time, center, Eigen::Matrix3d::Identity(), + return A.createFromArray(mass, pos, time, center, RowMatrix3d::Identity(), roundrobin, posvelrows); }, R"( @@ -1267,12 +1251,7 @@ void BasisFactoryClasses(py::module &m) py::arg("mass"), py::arg("pos"), py::arg("time"), py::arg("center") = Eigen::Vector3d::Zero(), py::arg("roundrobin") = true, py::arg("posvelrows") = true) - .def("initFromArray", - [](BasisClasses::BiorthBasis& A, - Eigen::Vector3d ctr, Eigen::Matrix3d rot) - { - return A.initFromArray(ctr, rot); - }, + .def("initFromArray", &BasisClasses::BiorthBasis::initFromArray, R"( Initialize coefficient accumulation @@ -1280,6 +1259,9 @@ void BasisFactoryClasses(py::module &m) ---------- center : list, default=[0, 0, 0] vector of center positions + rotation : numpy.ndarray + the rotation matrix for the basis functions + (default is identity matrix) Returns ------- @@ -1297,11 +1279,12 @@ void BasisFactoryClasses(py::module &m) phase-space sets that can not be stored in physical memory. )", py::arg("center") = Eigen::Vector3d::Zero(), - py::arg("rotation") = Eigen::Matrix3d::Identity()) + py::arg("rotation") = RowMatrix3d::Identity()) .def("addFromArray", - [](BasisClasses::BiorthBasis& A, Eigen::VectorXd& mass, RowMatrixXd& pos) + [](BasisClasses::BiorthBasis& A, + Eigen::VectorXd& mass, RowMatrixXd& pos) { - return A.addFromArray(mass, pos); + A.addFromArray(mass, pos); }, R"( Add particle contributions to coefficients @@ -2204,6 +2187,9 @@ void BasisFactoryClasses(py::module &m) the ParticleReader instance center : list, default=[0, 0, 0] an optional expansion center location + rotation : numpy.ndarray + the rotation matrix for the basis functions + (default is identity matrix) Returns ------- @@ -2212,13 +2198,8 @@ void BasisFactoryClasses(py::module &m) )", py::arg("reader"), py::arg("center") = Eigen::Vector3d::Zero(), - py::arg("rotation") = Eigen::Matrix3d::Identity()) - .def("initFromArray", - [](BasisClasses::FieldBasis& A, - Eigen::Vector3d center, Eigen::MatrixXd rotation) - { - return A.initFromArray(center, rotation); - }, + py::arg("rotation") = RowMatrix3d::Identity()) + .def("initFromArray", &BasisClasses::FieldBasis::initFromArray, R"( Initialize coefficient accumulation @@ -2245,11 +2226,11 @@ void BasisFactoryClasses(py::module &m) phase-space sets that can not be stored in physical memory. )", py::arg("center") = Eigen::Vector3d::Zero(), - py::arg("rotation") = Eigen::Matrix3d::Identity()) + py::arg("rotation") = RowMatrix3d::Identity()) .def("initFromArray", [](BasisClasses::FieldBasis& A, Eigen::Vector3d center) { - return A.initFromArray(center, Eigen::Matrix3d::Identity()); + return A.initFromArray(center, RowMatrix3d::Identity()); }, R"( Initialize coefficient accumulation @@ -2278,9 +2259,10 @@ void BasisFactoryClasses(py::module &m) )", py::arg("center") = Eigen::Vector3d::Zero()) .def("addFromArray", - [](BasisClasses::FieldBasis& A, Eigen::VectorXd& mass, RowMatrixXd& ps) + [](BasisClasses::FieldBasis& A, + Eigen::VectorXd& mass, RowMatrixXd& pos) { - return A.addFromArray(mass, ps); + return A.addFromArray(mass, pos); }, R"( Add particle contributions to coefficients @@ -2302,11 +2284,7 @@ void BasisFactoryClasses(py::module &m) makeFromArray : create coefficients contributions )", py::arg("mass"), py::arg("pos")) - .def("makeFromArray", - [](BasisClasses::FieldBasis& A, double time) - { - return A.makeFromArray(time); - }, + .def("makeFromArray", &BasisClasses::FieldBasis::makeFromArray, R"( Make the coefficients From bc0fe932b96792e3f2e25ee7fb4875ad0ff4b452 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 17 Aug 2025 19:35:31 -0400 Subject: [PATCH 15/18] Remove some debug output --- expui/BiorthBasis.cc | 3 --- 1 file changed, 3 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index b2b507717..093d946eb 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4016,9 +4016,6 @@ namespace BasisClasses Eigen::Map pp(p1.data(), 3), vv(v1.data(), 3); vv.setZero(); - if (myid==0) std::cout << "Center: " << coefctr << std::endl - << "Orient: " << coefrot << std::endl; - if (PosVelRows) { if (p.rows()<3) { std::ostringstream msg; From 188517e47ae51814798d19adf9711ec55acda8d6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 18 Aug 2025 10:51:07 -0400 Subject: [PATCH 16/18] Make rotation matrices row-major everywhere --- expui/CoefStruct.H | 10 ++++++---- pyEXP/CoefWrappers.cc | 2 +- src/Component.H | 6 ++++-- src/Cylinder.cc | 3 +++ src/Orient.H | 9 +++++---- src/PolarBasis.cc | 3 +++ src/SphericalBasis.cc | 3 +++ 7 files changed, 25 insertions(+), 11 deletions(-) diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index 88e1f2e4b..b7ee8ac7b 100644 --- a/expui/CoefStruct.H +++ b/expui/CoefStruct.H @@ -7,6 +7,8 @@ namespace CoefClasses { + using RowMatrix3d = Eigen::Matrix; + //! This is a general struct for any dimensionality class CoefStruct { @@ -27,10 +29,10 @@ namespace CoefClasses Eigen::VectorXcd store; //! Center data - Eigen::Vector3d ctr= Eigen::Vector3d::Zero(); + Eigen::Vector3d ctr = Eigen::Vector3d::Zero(); //! Orientation data - Eigen::Matrix3d rot= Eigen::Matrix3d::Identity(); + RowMatrix3d rot = RowMatrix3d::Identity(); //! Destructor virtual ~CoefStruct() {} @@ -72,10 +74,10 @@ namespace CoefClasses Eigen::Vector3d getCenter() { return ctr; } //! Set new rotation matrix - void setRotation(Eigen::Matrix3d& ROT) { rot = ROT; } + void setRotation(RowMatrix3d& ROT) { rot = ROT; } //! Read-only access to orientation - Eigen::Matrix3d getRotation() { return rot; } + RowMatrix3d getRotation() { return rot; } //! Set coefficient time (no copy) void setTime(double& STORE) diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index 48d344049..38e6beabe 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -811,7 +811,7 @@ void CoefficientClasses(py::module &m) { setCoefRotation : read-write access to the rotation matrix )") .def("setCoefRotation", - static_cast(&CoefStruct::setRotation), + static_cast(&CoefStruct::setRotation), py::arg("mat"), R"( Set the rotation matrix diff --git a/src/Component.H b/src/Component.H index 0edcea661..47e0ee84c 100644 --- a/src/Component.H +++ b/src/Component.H @@ -192,6 +192,8 @@ class Component friend class OrbTrace; friend class ComponentContainer; + using RowMatrix3d = Eigen::Matrix; + private: //! Parallel distribute and particle io @@ -767,9 +769,9 @@ public: } //! Get rotation matrix - Eigen::Matrix3d getRotation() + RowMatrix3d getRotation() { - Eigen::Matrix3d ret = Eigen::Matrix3d::Identity(); + RowMatrix3d ret = RowMatrix3d::Identity(); if ( (EJ & Orient::AXIS) and !EJdryrun) ret = orient->transformBody(); return ret; } diff --git a/src/Cylinder.cc b/src/Cylinder.cc index af8103e0d..b543c2ae4 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -1557,6 +1557,9 @@ void Cylinder::dump_coefs_h5(const std::string& file) // cur->ctr = component->getCenter(Component::Local | Component::Centered); + // Add the orientation + cur->rot = component->getRotation(); + // Check if file exists // if (std::filesystem::exists(file)) { diff --git a/src/Orient.H b/src/Orient.H index bb106bd51..4747344f3 100644 --- a/src/Orient.H +++ b/src/Orient.H @@ -88,7 +88,8 @@ Log file contents: */ class Orient { - typedef pair DV; + using DV = pair; + using RowMatrix3d = Eigen::Matrix; private: EL3 t; @@ -96,7 +97,7 @@ private: deque sumsA, sumsC; int keep, current; double damp; - Eigen::Matrix3d body, orig; + RowMatrix3d body, orig; Eigen::Vector3d axis, center, axis1, center1, center0, cenvel0; Eigen::Vector3d pseudo = Eigen::Vector3d::Zero(); Eigen::Vector3d omega = Eigen::Vector3d::Zero(); @@ -162,10 +163,10 @@ public: void accumulate(double time, Component* c); //! Return transformation to new body axes - Eigen::Matrix3d& transformBody(void) {return body;}; + RowMatrix3d& transformBody(void) {return body;}; //! Return transformation to original coordinates - Eigen::Matrix3d& transformOrig(void) {return orig;}; + RowMatrix3d& transformOrig(void) {return orig;}; //! Return current axis Eigen::Vector3d& currentAxis(void) {return axis;}; diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc index aae196aeb..e87b792fc 100644 --- a/src/PolarBasis.cc +++ b/src/PolarBasis.cc @@ -1822,6 +1822,9 @@ void PolarBasis::dump_coefs_h5(const std::string& file) // cur->ctr = component->getCenter(Component::Local | Component::Centered); + // Add the orientation + cur->rot = component->getRotation(); + // Check if file exists // if (std::filesystem::exists(file)) { diff --git a/src/SphericalBasis.cc b/src/SphericalBasis.cc index 4cde80982..74ff8feec 100644 --- a/src/SphericalBasis.cc +++ b/src/SphericalBasis.cc @@ -1906,6 +1906,9 @@ void SphericalBasis::dump_coefs_h5(const std::string& file) // cur->ctr = component->getCenter(Component::Local | Component::Centered); + // Add the orientation + cur->rot = component->getRotation(); + // Check if file exists // if (std::filesystem::exists(file)) { From 9a7d56d703f7f390f002c8e281cdafaf3081cfa1 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Mon, 29 Sep 2025 13:47:04 +0100 Subject: [PATCH 17/18] Update pyEXP/BasisWrappers.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index f2fd10a5a..6fede7399 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -298,8 +298,8 @@ void BasisFactoryClasses(py::module &m) } - virtual void initFromArray - (Eigen::Vector3d ctr, RowMatrix3d rot) { + void initFromArray + (Eigen::Vector3d ctr, RowMatrix3d rot) override { PYBIND11_OVERRIDE_PURE(void, Basis, initFromArray, ctr, rot); } From 8c48318a2605f82f49f8676068bd7590a10bcf69 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Mon, 29 Sep 2025 09:45:48 -0400 Subject: [PATCH 18/18] Update pyEXP/BasisWrappers.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 6fede7399..9f272135b 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2262,7 +2262,7 @@ void BasisFactoryClasses(py::module &m) [](BasisClasses::FieldBasis& A, Eigen::VectorXd& mass, RowMatrixXd& pos) { - return A.addFromArray(mass, pos); + A.addFromArray(mass, pos); }, R"( Add particle contributions to coefficients