rp/solvers/Mixture.h

#ifndef RP_MIXTURE_H
#define RP_MIXTURE_H
 
/*
* PSCF - Polymer Self-Consistent Field
*
* Copyright 2015 - 2025, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include <pscf/solvers/MixtureTmpl.h>   // base class template
#include <util/containers/FSArray.h>    // member template (stress_)
#include <iostream>
 
// Forward declarations
namespace Util {
   template <typename T> class DArray;
}
namespace Pscf {
   template <int D> class Mesh;
   namespace Prdc {
      template <int D> class UnitCell;
   }
}
 
namespace Pscf {
namespace Rp {
 
   using namespace Util;
   using namespace Prdc;

   template <int D, class T>
   class Mixture : public 
      MixtureTmpl<typename T::Polymer, typename T::Solvent, double>
   {
 
   public:
 
      // Public type name aliases

      using PolymerT = typename T::Polymer;

      using SolventT = typename T::Solvent;

      using BlockT = typename T::Block;

      using PropagatorT = typename T::Propagator;

      using FieldT = typename T::RField;

      using FFTT = typename T::FFT;

      using WaveListT = typename T::WaveList;

      using FieldIoT = typename T::FieldIo;

      using MixtureTmplT = MixtureTmpl<PolymerT, SolventT, double>;

      using MixtureBaseT = typename MixtureTmplT::MixtureBaseT;
 
      // Public member functions


      virtual void readParameters(std::istream& in);

      void associate(Mesh<D> const & mesh,
                     FFTT const & fft,
                     UnitCell<D> const & cell,
                     WaveListT& waveList);

      void setFieldIo(FieldIoT const & fieldIo);

      void allocate();


      void compute(DArray<FieldT> const & wFields,
                   DArray<FieldT>& cFields,
                   double phiTot = 1.0);

      void setIsSymmetric(bool isSymmetric);

      void computeStress(double phiTot = 1.0);

      bool hasStress() const;

      double stress(int parameterId) const;


      void setKuhn(int monomerId, double kuhn);

      void clearUnitCellData();


      void createBlockCRGrid(DArray<FieldT>& blockCFields) const;

      void writeBlockCRGrid(std::string const & filename) const;

      
      void writeQSlice(std::string const & filename,
                       int polymerId, int blockId,
                       int directionId, int segmentId)  const;

      void writeQTail(std::string const & filename, int polymerId,
                      int blockId, int directionId)  const;

      void writeQ(std::string const & filename, int polymerId,
                  int blockId, int directionId)  const;

      void writeQAll(std::string const & basename);


      void writeStress(std::ostream& out) const;

 
      // Inherited non-dependent public member functions
      using MixtureTmplT::polymer;
      using MixtureTmplT::polymerSpecies;
      using MixtureTmplT::solvent;
      using MixtureTmplT::solventSpecies;
      using MixtureBaseT::nMonomer;
      using MixtureBaseT::monomer;
      using MixtureBaseT::nPolymer;
      using MixtureBaseT::nSolvent;
      using MixtureBaseT::nBlock;
      using MixtureBaseT::vMonomer;
      using MixtureBaseT::isCanonical;
 
   protected:

      Mixture();

      ~Mixture();

      Mesh<D> const & mesh() const
      {  return *meshPtr_; }

      UnitCell<D> const & unitCell() const
      {  return *unitCellPtr_; }

      FieldIoT const & fieldIo() const
      {  return *fieldIoPtr_; }

      double ds() const
      {  return ds_; }
 
   private:
 
      // Private member data

      FSArray<double, 6> stress_;

      double ds_;

      Mesh<D> const * meshPtr_;
 
      // Pointer to associated UnitCell<D> object 
      UnitCell<D> const * unitCellPtr_;
 
      // Pointer to associated FieldIoT object 
      FieldIoT const * fieldIoPtr_;

      int nParam_;

      bool hasStress_;
 
      // Set true iff the w fields used in the MDE are symmetric
      bool isSymmetric_;
 
      // Private member functions (pure virtual)

      virtual void allocateBlocks() = 0; 
 
   };
 
   // Public inline member functions
 
   /*
   * Has the stress been computed for the current w fields?
   */
   template <int D, class T>
   inline bool Mixture<D,T>::hasStress() const
   {  return hasStress_; }
 
   /*
   * Get derivative of free energy w/ respect to a unit cell parameter.
   */
   template <int D, class T>
   inline double Mixture<D,T>::stress(int parameterId) const
   {
      UTIL_CHECK(hasStress_);
      UTIL_CHECK(parameterId < nParam_);
      return stress_[parameterId];
   }
 
} // namespace Rp
} // namespace Pscf
#endif