MixtureReal.h

#ifndef PRDC_MIXTURE_REAL_H
#define PRDC_MIXTURE_REAL_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 Prdc {
 
   using namespace Util;

   template <int D, class PT, class ST>
   class MixtureReal : public MixtureTmpl<PT, ST>
   {
 
   public:
 
      // Public type name aliases

      using MixtureTmplT = MixtureTmpl<PT,ST>;

      using typename MixtureTmplT::SolventT;

      using typename MixtureTmplT::PolymerT;

      using typename MixtureTmplT::BlockT;

      using typename MixtureTmplT::PropagatorT;

      using FieldT = typename PropagatorT::FieldT;

      using FFTT = typename BlockT::FFTT;

      using WaveListT = typename BlockT::WaveListT;

      using FieldIoT = typename BlockT::FieldIoT;
 
      // Public member functions


      MixtureReal();

      ~MixtureReal();

      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 public member functions
      using MixtureTmplT::polymer;
      using MixtureTmplT::polymerSpecies;
      using MixtureTmplT::solvent;
      using MixtureTmplT::solventSpecies;
      using MixtureBase::nMonomer;
      using MixtureBase::monomer;
      using MixtureBase::nPolymer;
      using MixtureBase::nSolvent;
      using MixtureBase::nBlock;
      using MixtureBase::vMonomer;
      using MixtureBase::isCanonical;
 
   protected:

      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_; }
 
      // Inherited protected member functions
      using ParamComposite::setClassName;
      using ParamComposite::read;
      using ParamComposite::readOptional;
 
   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

      virtual void eqS(FieldT& A, double s) const = 0;

      virtual void addEqV(FieldT& A, FieldT const & B) const = 0;

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