PolymerTmpl.h

#ifndef PSCF_POLYMER_TMPL_H
#define PSCF_POLYMER_TMPL_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 <pscf/chem/PolymerSpecies.h>    // base class
#include <util/containers/DArray.h>      // member template
 
// Classes used in implementation
#include <pscf/chem/Vertex.h>           
#include <pscf/chem/PolymerType.h>       
#include <util/containers/Pair.h>       
 
#include <cmath>
 
namespace Pscf
{
 
   class Edge;
   using namespace Util;
 
   template <class Block>
   class PolymerTmpl : public PolymerSpecies
   {
 
   public:
 
      // Modified diffusion equation solver for one block.
      typedef typename Block::Propagator Propagator;
 
      // Monomer concentration field.
      typedef typename Propagator::CField CField;
 
      // Chemical potential field.
      typedef typename Propagator::WField WField;
 
      PolymerTmpl();
 
      ~PolymerTmpl();
 
      virtual void readParameters(std::istream& in);
 
      virtual void solve(double phiTot = 1.0);
 
 
      virtual Edge& edge(int id);
 
      virtual Edge const& edge(int id) const;
 
      Block& block(int id);
 
      Block const& block(int id) const ;
 
      Propagator& propagator(int blockId, int directionId);
 
      Propagator const & propagator(int blockId, int directionId) const;
 
      Propagator& propagator(int id);
 
 
      // Inherited public members
      using PolymerSpecies::edge;
      using PolymerSpecies::vertex;
      using PolymerSpecies::propagatorId;
      using PolymerSpecies::path;
      using PolymerSpecies::nBlock;
      using PolymerSpecies::nVertex;
      using PolymerSpecies::nPropagator;
      using PolymerSpecies::length;
      using PolymerSpecies::type;
 
      using Species::phi;
      using Species::mu;
      using Species::q;
      using Species::ensemble;
 
   protected:
 
      virtual void allocateBlocks();
 
      virtual void readBlocks(std::istream& in);
 
   private:
 
      DArray<Block> blocks_;
 
   };
 
   /*
   * Get a specified Edge (block descriptor) by non-const reference.
   */
   template <class Block>
   inline Edge& PolymerTmpl<Block>::edge(int id)
   {  return blocks_[id]; }
 
   /*
   * Get a specified Edge (block descriptor) by const reference.
   */
   template <class Block>
   inline 
   Edge const & PolymerTmpl<Block>::edge(int id) const
   {  return blocks_[id]; }
 
   /*
   * Get a specified Block solver by non-const reference.
   */
   template <class Block>
   inline 
   Block& PolymerTmpl<Block>::block(int id)
   {  return blocks_[id]; }
 
   /*
   * Get a specified Block solver by const reference.
   */
   template <class Block>
   inline 
   Block const & PolymerTmpl<Block>::block(int id) const
   {  return blocks_[id]; }
 
   /*
   * Get a Propagator indexed by block and direction (non-const).
   */
   template <class Block>
   inline
   typename Block::Propagator&
   PolymerTmpl<Block>::propagator(int blockId, int directionId)
   {  return block(blockId).propagator(directionId); }
 
   /*
   * Get a Propagator indexed by block and direction (const).
   */
   template <class Block>
   inline
   typename Block::Propagator const &
   PolymerTmpl<Block>::propagator(int blockId, int directionId) const
   {  return block(blockId).propagator(directionId); }
 
   /*
   * Get a propagator indexed in order of computation.
   */
   template <class Block>
   inline
   typename Block::Propagator& PolymerTmpl<Block>::propagator(int id)
   {
      Pair<int> propId = propagatorId(id);
      return propagator(propId[0], propId[1]);
   }
 
   // Non-inline functions
 
   /*
   * Constructor.
   */
   template <class Block>
   PolymerTmpl<Block>::PolymerTmpl()
    : PolymerSpecies(),
      blocks_() 
   {  setClassName("PolymerTmpl"); }
 
   /*
   * Destructor.
   */
   template <class Block>
   PolymerTmpl<Block>::~PolymerTmpl()
   {}
 
   template <class Block>
   void PolymerTmpl<Block>::allocateBlocks()
   {  blocks_.allocate(nBlock()); }
 
   template <class Block>
   void PolymerTmpl<Block>::readBlocks(std::istream& in)
   {  readDArray<Block>(in, "blocks", blocks_, nBlock()); }
 
   template <class Block>
   void PolymerTmpl<Block>::readParameters(std::istream& in)
   {
 
      PolymerSpecies::readParameters(in);
 
      // Set sources for all propagators
      Vertex const * vertexPtr = 0;
      Propagator const * sourcePtr = 0;
      Propagator * propagatorPtr = 0;
      Pair<int> propagatorId;
      int blockId, directionId, vertexId, i;
      for (blockId = 0; blockId < nBlock(); ++blockId) {
         // Add sources
         for (directionId = 0; directionId < 2; ++directionId) {
            vertexId = block(blockId).vertexId(directionId);
            vertexPtr = &vertex(vertexId);
            propagatorPtr = &block(blockId).propagator(directionId);
            for (i = 0; i < vertexPtr->size(); ++i) {
               propagatorId = vertexPtr->inPropagatorId(i);
               if (propagatorId[0] == blockId) {
                  UTIL_CHECK(propagatorId[1] != directionId);
               } else {
                  sourcePtr =
                     &block(propagatorId[0]).propagator(propagatorId[1]);
                  propagatorPtr->addSource(*sourcePtr);
               }
            }
         }
 
      }
 
   }
 
   /*
   * Compute a solution to the MDE and block concentrations.
   */
   template <class Block>
   void PolymerTmpl<Block>::solve(double phiTot)
   {
 
      // Clear all propagators
      for (int j = 0; j < nPropagator(); ++j) {
         propagator(j).setIsSolved(false);
      }
 
      // Solve modified diffusion equation for all propagators in
      // the order specified by function makePlan.
      for (int j = 0; j < nPropagator(); ++j) {
         UTIL_CHECK(propagator(j).isReady());
         propagator(j).solve();
      }
 
      // Compute molecular partition function q_
      q_ = block(0).propagator(0).computeQ();
 
      // The Propagator::computeQ function returns a spatial average.
      // Correct for partial occupation of the unit cell.
      q_ = q_/phiTot;
 
      // Compute mu_ or phi_, depending on ensemble
      if (ensemble() == Species::Closed) {
         mu_ = log(phi_/q_);
      } else
      if (ensemble() == Species::Open) {
         phi_ = exp(mu_)*q_;
      }
 
      // Compute block concentration fields
      double prefactor = phi_ / ( q_ * length() );
      for (int i = 0; i < nBlock(); ++i) {
         block(i).computeConcentration(prefactor);
      }
 
   }
 
}
#endif