.

include/sctl/fft_wrapper.hpp

@@ -104,8 +104,17 @@ template <class ValueType, class FFT_Derived> class FFT_Generic {
     dim[0] = 0;
     dim[1] = 0;
   }
-  FFT_Generic(const FFT_Generic&) = delete;
-  FFT_Generic& operator=(const FFT_Generic&) = delete;
+
+  FFT_Generic(const FFT_Generic&) {
+    dim[0]=0;
+    dim[1]=0;
+  }
+
+  FFT_Generic& operator=(const FFT_Generic&) {
+    dim[0]=0;
+    dim[1]=0;
+    return *this;
+  };
 
   Long Dim(Integer i) const {
     return dim[i];

include/sctl/sph_harm.hpp

@@ -63,6 +63,8 @@ template <class Real> class SphericalHarmonics{
       SphericalHarmonics<Real>::WriteVTK("test", nullptr, &Xcoeff, sctl::SHCArrange::ROW_MAJOR, p, 32);
     }
 
+    static void Clear() { MatrixStore().Resize(0); }
+
   private:
 
     // Probably don't work anymore, need to be updated :(
@@ -106,6 +108,10 @@ template <class Real> class SphericalHarmonics{
     struct MatrixStorage{
       MatrixStorage(){
         const Long size = SCTL_SHMAXDEG;
+        Resize(size);
+      }
+
+      void Resize(Long size){
         Qx_ .resize(size);
         Qw_ .resize(size);
         Sw_ .resize(size);
@@ -116,7 +122,11 @@ template <class Real> class SphericalHarmonics{
         Mr_ .resize(size);
         Mfinv_ .resize(size*size);
         Mlinv_ .resize(size*size);
+
+        Mfft_.resize(size);
+        Mfftinv_.resize(size);
       }
+
       std::vector<Vector<Real>> Qx_;
       std::vector<Vector<Real>> Qw_;
       std::vector<Vector<Real>> Sw_;
@@ -128,8 +138,8 @@ template <class Real> class SphericalHarmonics{
       std::vector<Matrix<Real>> Mfinv_ ;
       std::vector<std::vector<Matrix<Real>>> Mlinv_ ;
 
-      StaticArray<FFT<Real>, SCTL_SHMAXDEG> Mfft_ ;
-      StaticArray<FFT<Real>, SCTL_SHMAXDEG> Mfftinv_ ;
+      std::vector<FFT<Real>> Mfft_;
+      std::vector<FFT<Real>> Mfftinv_;
     };
     static MatrixStorage& MatrixStore(){
       static MatrixStorage storage;
@@ -137,6 +147,8 @@ template <class Real> class SphericalHarmonics{
     }
 };
 
+template class SphericalHarmonics<double>;
+
 }  // end namespace
 
 #include SCTL_INCLUDE(sph_harm.txx)

include/sctl/sph_harm.txx

@@ -885,8 +885,11 @@ template <class Real> const FFT<Real>& SphericalHarmonics<Real>::OpFourier(Long
   assert(Np<SCTL_SHMAXDEG);
   auto& Mf =MatrixStore().Mfftinv_ [Np];
   if(!Mf.Dim(0)){
-    StaticArray<Long,1> fft_dim = {Np};
-    Mf.Setup(FFT_Type::C2R, 1, Vector<Long>(1,fft_dim,false));
+    #pragma omp critical SCTL_FFT_PLAN0
+    if(!Mf.Dim(0)){
+      StaticArray<Long,1> fft_dim = {Np};
+      Mf.Setup(FFT_Type::C2R, 1, Vector<Long>(1,fft_dim,false));
+    }
   }
   return Mf;
 }
@@ -895,8 +898,11 @@ template <class Real> const FFT<Real>& SphericalHarmonics<Real>::OpFourierInv(Lo
   assert(Np<SCTL_SHMAXDEG);
   auto& Mf =MatrixStore().Mfft_ [Np];
   if(!Mf.Dim(0)){
-    StaticArray<Long,1> fft_dim = {Np};
-    Mf.Setup(FFT_Type::R2C, 1, Vector<Long>(1,fft_dim,false));
+    #pragma omp critical SCTL_FFT_PLAN1
+    if(!Mf.Dim(0)){
+      StaticArray<Long,1> fft_dim = {Np};
+      Mf.Setup(FFT_Type::R2C, 1, Vector<Long>(1,fft_dim,false));
+    }
   }
   return Mf;
 }