#ifndef _SCTL_SPH_HARM_HPP_ #define _SCTL_SPH_HARM_HPP_ #define SCTL_SHMAXDEG 1024 #include SCTL_INCLUDE(matrix.hpp) #include SCTL_INCLUDE(fft_wrapper.hpp) #include SCTL_INCLUDE(common.hpp) namespace SCTL_NAMESPACE { enum class SHCArrange { // (p+1) x (p+1) complex elements in row-major order. // A : { A(0,0), A(0,1), ... A(0,p), A(1,0), ... A(p,p) } // where, A(n,m) = { Ar(n,m), Ai(n,m) } (real and imaginary parts) ALL, // (p+1)(p+2)/2 complex elements in row-major order (lower triangular part) // A : { A(0,0), A(1,0), A(1,1), A(2,0), A(2,1), A(2,2), ... A(p,p) } // where, A(n,m) = { Ar(n,m), Ai(n,m) } (real and imaginary parts) ROW_MAJOR, // (p+1)(p+1) real elements in col-major order (non-zero lower triangular part) // A : { Ar(0,0), Ar(1,0), ... Ar(p,0), Ar(1,1), ... Ar(p,1), Ai(1,1), ... Ai(p,1), ..., Ar(p,p), Ai(p,p) // where, A(n,m) = { Ar(n,m), Ai(n,m) } (real and imaginary parts) COL_MAJOR_NONZERO }; template class SphericalHarmonics{ static constexpr Integer COORD_DIM = 3; public: // Scalar Spherical Harmonics static void Grid2SHC(const Vector& X_in, Long Nt_in, Long Np_in, Long p_out, Vector& S_out, SHCArrange arrange_out); static void SHC2Grid(const Vector& S_in, SHCArrange arrange_in, Long p_in, Long Nt_out, Long Np_out, Vector* X_out, Vector* X_theta_out=nullptr, Vector* X_phi_out=nullptr); static void SHCEval(const Vector& S_in, SHCArrange arrange_in, Long p_in, const Vector& theta_phi_in, Vector& X_out); static void SHC2Pole(const Vector& S_in, SHCArrange arrange_in, Long p_in, Vector& P_out); static void WriteVTK(const char* fname, const Vector* S, const Vector* f_val, SHCArrange arrange, Long p_in, Long p_out, Real period=0, const Comm& comm = Comm::World()); // Vector Spherical Harmonics static void Grid2VecSHC(const Vector& X_in, Long Nt_in, Long Np_in, Long p_out, Vector& S_out, SHCArrange arrange_out); static void VecSHC2Grid(const Vector& S_in, SHCArrange arrange_in, Long p_in, Long Nt_out, Long Np_out, Vector* X_out); static void VecSHCEval(const Vector& S_in, SHCArrange arrange_in, Long p_in, const Vector& theta_phi_in, Vector& X_out); static void StokesEvalSL(const Vector& S_in, SHCArrange arrange_in, Long p_in, const Vector& theta_phi_in, Vector& X_out); static void test() { int p = 4; int dof = 6; int Nt = p+1, Np = 2*p+1; auto print_coeff = [&](Vector S) { Long idx=0; for (Long k=0;k(2*n+2, S.begin()+idx); idx+=2*n+2; } } std::cout<<'\n'; }; Vector r_theta_phi, theta_phi; { // Set r_theta_phi, theta_phi Vector leg_nodes = LegendreNodes(Nt-1); for (Long i=0;i() / Np); theta_phi.PushBack(leg_nodes[i]); theta_phi.PushBack(j * 2 * const_pi() / Np); } } } int Ncoeff = (p + 1) * (p + 1); Vector Xcoeff(dof * Ncoeff), Xgrid; for (int i=0;i(Nt*dof, Np, Xgrid.begin())<<'\n'; { Vector val; VecSHCEval(Xcoeff, sctl::SHCArrange::COL_MAJOR_NONZERO, p, theta_phi, val); //StokesEvalSL(Xcoeff, sctl::SHCArrange::COL_MAJOR_NONZERO, p, r_theta_phi, val); Matrix(dof, val.Dim()/dof, val.begin(), false) = Matrix(val.Dim()/dof, dof, val.begin()).Transpose(); std::cout<(val.Dim()/Np, Np, val.begin()) - Matrix(Nt*dof, Np, Xgrid.begin())<<'\n'; } Grid2VecSHC(Xgrid, Nt, Np, p, Xcoeff, sctl::SHCArrange::ROW_MAJOR); print_coeff(Xcoeff); //SphericalHarmonics::WriteVTK("test", nullptr, &Xcoeff, sctl::SHCArrange::ROW_MAJOR, p, 32); Clear(); } static void Clear() { MatrixStore().Resize(0); } private: // Probably don't work anymore, need to be updated :( static void SHC2GridTranspose(const Vector& X, Long p0, Long p1, Vector& S); static void RotateAll(const Vector& S, Long p0, Long dof, Vector& S_); static void RotateTranspose(const Vector& S_, Long p0, Long dof, Vector& S); static void StokesSingularInteg(const Vector& S, Long p0, Long p1, Vector* SLMatrix=nullptr, Vector* DLMatrix=nullptr); static void Grid2SHC_(const Vector& X, Long Nt, Long Np, Long p, Vector& B1); static void SHCArrange0(const Vector& B1, Long p, Vector& S, SHCArrange arrange); static void SHC2Grid_(const Vector& S, Long p, Long Nt, Long Np, Vector* X, Vector* X_theta=nullptr, Vector* X_phi=nullptr); static void SHCArrange1(const Vector& S_in, SHCArrange arrange_out, Long p, Vector& S_out); /** * \brief Computes all the Associated Legendre Polynomials (normalized) up to the specified degree. * \param[in] degree The degree up to which the Legendre polynomials have to be computed. * \param[in] X The input values for which the polynomials have to be computed. * \param[in] N The number of input points. * \param[out] poly_val The output array of size (degree+1)*(degree+2)*N/2 containing the computed polynomial values. * The output values are in the order: * P(n,m)[i] => {P(0,0)[0], P(0,0)[1], ..., P(0,0)[N-1], P(1,0)[0], ..., P(1,0)[N-1], * P(2,0)[0], ..., P(degree,0)[N-1], P(1,1)[0], ...,P(2,1)[0], ..., P(degree,degree)[N-1]} */ static void LegPoly(Vector& poly_val, const Vector& X, Long degree); static void LegPolyDeriv(Vector& poly_val, const Vector& X, Long degree); static const Vector& LegendreNodes(Long p1); static const Vector& LegendreWeights(Long p1); static const Vector& SingularWeights(Long p1); static const Matrix& MatFourier(Long p0, Long p1); static const Matrix& MatFourierInv(Long p0, Long p1); static const Matrix& MatFourierGrad(Long p0, Long p1); static const FFT& OpFourier(Long Np); static const FFT& OpFourierInv(Long Np); static const std::vector>& MatLegendre(Long p0, Long p1); static const std::vector>& MatLegendreInv(Long p0, Long p1); static const std::vector>& MatLegendreGrad(Long p0, Long p1); // Evaluate all Spherical Harmonic basis functions up to order p at (theta, phi) coordinates. static void SHBasisEval(Long p, const Vector& cos_theta_phi, Matrix& M); static void VecSHBasisEval(Long p, const Vector& cos_theta_phi, Matrix& M); static const std::vector>& MatRotate(Long p0); template static void StokesSingularInteg_(const Vector& X0, Long p0, Long p1, Vector& SL, Vector& DL); struct MatrixStorage{ MatrixStorage(){ const Long size = SCTL_SHMAXDEG; Resize(size); } void Resize(Long size){ Qx_ .resize(size); Qw_ .resize(size); Sw_ .resize(size); Mf_ .resize(size*size); Mdf_.resize(size*size); Ml_ .resize(size*size); Mdl_.resize(size*size); Mr_ .resize(size); Mfinv_ .resize(size*size); Mlinv_ .resize(size*size); Mfft_.resize(size); Mfftinv_.resize(size); } std::vector> Qx_; std::vector> Qw_; std::vector> Sw_; std::vector> Mf_ ; std::vector> Mdf_; std::vector>> Ml_ ; std::vector>> Mdl_; std::vector>> Mr_; std::vector> Mfinv_ ; std::vector>> Mlinv_ ; std::vector> Mfft_; std::vector> Mfftinv_; }; static MatrixStorage& MatrixStore(){ static MatrixStorage storage; return storage; } }; template class SphericalHarmonics; } // end namespace #include SCTL_INCLUDE(sph_harm.txx) #endif // _SCTL_SPH_HARM_HPP_