rpc/solvers/Block.h

#ifndef RPC_BLOCK_H
#define RPC_BLOCK_H
 
/*
* PSCF - Polymer Self-Consistent Field Theory
*
* Copyright 2016 - 2022, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include "Propagator.h"                   // base class argument
#include <prdc/cpu/FFT.h>                 // member
#include <prdc/cpu/RField.h>              // member
#include <prdc/cpu/RFieldDft.h>           // member
#include <prdc/crystal/UnitCell.h>        // member
#include <pscf/solvers/BlockTmpl.h>       // base class template
#include <pscf/mesh/Mesh.h>               // member
#include <util/containers/FArray.h>       // member template
#include <util/containers/DMatrix.h>      // member template
 
namespace Pscf {
   template <int D> class Mesh;
   namespace Prdc{
     template <int D> class UnitCell;
   }
}
 
namespace Pscf {
namespace Rpc {
 
   using namespace Util;
   using namespace Pscf::Prdc;
   using namespace Pscf::Prdc::Cpu;
 
   template <int D>
   class Block : public BlockTmpl< Propagator<D> >
   {
 
   public:
 
      Block();
 
      ~Block();
 
      void associate(Mesh<D> const& mesh, 
                     FFT<D> const& fft, 
                     UnitCell<D> const& cell);
 
      void allocate(double ds);
 
      void clearUnitCellData();
 
      void setLength(double newLength);
 
      void setKuhn(double kuhn);
 
      void setupSolver(RField<D> const & w);
 
      void step(RField<D> const & q, RField<D>& qNew);
 
      void computeConcentration(double prefactor);
 
      void computeStress(double prefactor);
 
      Mesh<D> const & mesh() const;
 
      FFT<D> const & fft() const;
 
      double ds() const;
 
      int ns() const;
 
      double stress(int n) const;
 
      // Functions with non-dependent names from BlockTmpl< Propagator<D> >
      using BlockTmpl< Propagator<D> >::setKuhn;
      using BlockTmpl< Propagator<D> >::propagator;
      using BlockTmpl< Propagator<D> >::cField;
      using BlockTmpl< Propagator<D> >::length;
      using BlockTmpl< Propagator<D> >::kuhn;
 
      // Functions with non-dependent names from Edge
      using Edge::setId;
      using Edge::setVertexIds;
      using Edge::setMonomerId;
      using Edge::setLength;
      using Edge::id;
      using Edge::monomerId;
      using Edge::vertexIds;
      using Edge::vertexId;
      using Edge::length;
 
   private:
 
      // Matrix to store derivatives of plane waves
      DMatrix<double> dGsq_;
 
      // Stress arising from this block
      FSArray<double, 6> stress_;
 
      // Fourier transform plan
      // FFT<D> fft_;
 
      // Array of elements containing exp(-K^2 b^2 ds/6)
      RField<D> expKsq_;
 
      // Array of elements containing exp(-W[i] ds/2)
      RField<D> expW_;
 
      // Array of elements containing exp(-K^2 b^2 ds/(6*2))
      RField<D> expKsq2_;
 
      // Array of elements containing exp(-W[i] (ds/2)*0.5)
      RField<D> expW2_;
 
      // Work array for real-space field (step size ds)
      RField<D> qr_;
 
      // Work array for real-space field (step size ds/2)
      RField<D> qr2_;
 
      // Work array for wavevector space field (step size ds)
      RFieldDft<D> qk_;
 
      // Work array for wavevector space field (step size ds/2)
      RFieldDft<D> qk2_;
 
      // Pointer to associated Mesh<D> object
      Mesh<D> const* meshPtr_;
 
      // Pointer to associated FFT<D> object
      FFT<D> const* fftPtr_;
 
      // Pointer to associated UnitCell<D> object
      UnitCell<D> const* unitCellPtr_;
 
      // Dimensions of wavevector mesh in real-to-complex transform
      IntVec<D> kMeshDimensions_;
 
      // Number of wavevectors in wavevector mesh 
      int kSize_;
 
      // Contour length step size (actual step size for this block)
      double ds_;
 
      // Contour length step size (value input in param file)
      double dsTarget_;
 
      // Number of contour grid points = # of contour steps + 1
      int ns_;
 
      // Have arrays been allocated ?
      bool isAllocated_;
 
      // Are expKsq_ arrays up to date ? (initialize false)
      bool hasExpKsq_;
 
      UnitCell<D> const & unitCell() const
      {  return *unitCellPtr_; }
 
      void computedGsq();
 
      void computeExpKsq();
 
   };
 
   // Inline member functions
 
   template <int D>
   inline int Block<D>::ns() const
   {  return ns_; }
 
   template <int D>
   inline double Block<D>::ds() const
   {  return ds_; }
 
   template <int D>
   inline double Block<D>::stress(int n) const
   {  return stress_[n]; }
 
   template <int D>
   inline Mesh<D> const & Block<D>::mesh() const
   {
      UTIL_ASSERT(meshPtr_);
      return *meshPtr_;
   }
 
   template <int D>
   inline FFT<D> const & Block<D>::fft() const
   {
      UTIL_ASSERT(fftPtr_);
      return *fftPtr_;
      // return fft_;
   }
 
   #ifndef RPC_BLOCK_TPP
   // Suppresse implicit instantiation
   extern template class Block<1>;
   extern template class Block<2>;
   extern template class Block<3>;
   #endif
 
}
}
#endif