SpaceSymmetry.h

#ifndef PRDC_SPACE_SYMMETRY_H
#define PRDC_SPACE_SYMMETRY_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/math/IntVec.h>
#include <util/math/Rational.h>
#include <util/containers/FMatrix.h>
#include <util/containers/FArray.h>
#include <util/format/Int.h>
 
#include <iostream>
 
namespace Pscf {
namespace Prdc {
 
   using namespace Util;
 
   template <int D> class SpaceSymmetry;
 
   // Non-member function declarations

   template <int D>
   bool operator == (const SpaceSymmetry<D>& A, const SpaceSymmetry<D>& B);

   template <int D>
   bool operator != (const SpaceSymmetry<D>& A, const SpaceSymmetry<D>& B);

   template <int D>
   SpaceSymmetry<D> 
   operator * (const SpaceSymmetry<D>& A, const SpaceSymmetry<D>& B);

   template <int D>
   IntVec<D> operator * (const SpaceSymmetry<D>& S, const IntVec<D>& V);

   template <int D>
   IntVec<D> operator * (const IntVec<D>& V, const SpaceSymmetry<D>& S);
 
   template <int D>
   std::ostream& operator << (std::ostream& out, const SpaceSymmetry<D>& A);
 
   template <int D>
   std::istream& operator >> (std::istream& in, SpaceSymmetry<D>& A);

   template <int D>
   class SpaceSymmetry 
   {
 
   public:

      typedef FMatrix<int, D, D> Rotation;

      typedef FArray<Rational, D> Translation;

      SpaceSymmetry();

      SpaceSymmetry(const SpaceSymmetry<D>& other);

      SpaceSymmetry<D>& operator = (const SpaceSymmetry<D>& other);

      void normalize();

      void shiftOrigin(Translation const & origin);

      SpaceSymmetry<D> inverse() const;

      SpaceSymmetry<D>::Rotation inverseRotation() const;

      int determinant() const;

      int& R(int i, int j);
 
      int R(int i, int j) const;

      Rational& t(int i);
 
      Rational t(int i) const;
 
      // Static method

      static const SpaceSymmetry<D>& identity();
 
   private:

      Rotation R_;

      Translation t_;
 
      // Static member variables

      static SpaceSymmetry<D> identity_;

      static bool hasIdentity_;

      static void makeIdentity();
 
   // friends:
 
      friend 
      bool operator == <> (const SpaceSymmetry<D>& A, 
                           const SpaceSymmetry<D>& B);
 
      friend 
      bool operator != <> (const SpaceSymmetry<D>& A, 
                           const SpaceSymmetry<D>& B);
 
      friend 
      SpaceSymmetry<D>
      operator * <> (const SpaceSymmetry<D>& A, 
                     const SpaceSymmetry<D>& B);
 
      friend 
      IntVec<D> operator * <> (const IntVec<D>& V, 
                               const SpaceSymmetry<D>& S);
 
      friend 
      IntVec<D> operator * <>(const SpaceSymmetry<D>& S, const IntVec<D>& V);
 
      friend 
      std::ostream& operator << <> (std::ostream& out, 
                                    const SpaceSymmetry<D>& A);
 
      friend 
      std::istream& operator >> <> (std::istream& in, 
                                    SpaceSymmetry<D>& A);
 
   };
 
   // Explicit specialization of members
 
   template <> 
   SpaceSymmetry<1>::Rotation SpaceSymmetry<1>::inverseRotation() const;
 
   template <> 
   SpaceSymmetry<2>::Rotation SpaceSymmetry<2>::inverseRotation() const;
 
   template <> 
   SpaceSymmetry<3>::Rotation SpaceSymmetry<3>::inverseRotation() const;
 
   template <> 
   int SpaceSymmetry<1>::determinant() const;
 
   template <> 
   int SpaceSymmetry<2>::determinant() const;
 
   template <> 
   int SpaceSymmetry<3>::determinant() const;
 
   // Inline method definitions
 
   template <int D>
   inline 
   bool operator != (const SpaceSymmetry<D>& A, const SpaceSymmetry<D>& B)
   {  return !(A == B); }
 
   /*
   * Return an element of the matrix by reference.
   */
   template <int D>
   inline 
   int& SpaceSymmetry<D>::R(int i, int j)
   {  return R_(i, j); }
 
   /*
   * Return an element of the matrix by value
   */
   template <int D>
   inline 
   int SpaceSymmetry<D>::R(int i, int j) const
   {  return R_(i, j); }
 
   /*
   * Return an element of the translation vector by reference.
   */
   template <int D>
   inline 
   Rational& SpaceSymmetry<D>::t(int i)
   {  return t_[i]; }
 
   /*
   * Return an element of the translation vector by value
   */
   template <int D>
   inline 
   Rational SpaceSymmetry<D>::t(int i) const
   {  return t_[i]; }
 
   /*
   * Return the identity symmetry operation.
   */
   template <int D>
   inline 
   const SpaceSymmetry<D>& SpaceSymmetry<D>::identity()
   {
      if (!hasIdentity_) makeIdentity();
      return identity_;
   }
 
   // Friend function template definitions
 
   /*
   * Equality operator.
   */
   template <int D>
   bool operator == (const SpaceSymmetry<D>& A, const SpaceSymmetry<D>& B)
   {
      int i, j;
      for (i = 0; i < D; ++i) {
         if (A.t_[i] != B.t_[i]) {
            return false;
         }
         for (j = 0; j < D; ++j) {
            if (A.R_(i, j) != B.R_(i,j)) {
               return false;
            }
         }
      }
      return true;
   }
 
   /*
   * Group multipication operator for SpaceSymmetry<D> objects.
   */
   template <int D>
   SpaceSymmetry<D> 
   operator * (const SpaceSymmetry<D>& A, const SpaceSymmetry<D>& B)
   {
      SpaceSymmetry<D> C;
      int i, j, k;
 
      // Compute rotation matrix (matrix product)
      for (i = 0; i < D; ++i) {
         for (j = 0; j < D; ++j) {
            C.R_(i, j) = 0;
            for (k = 0; k < D; ++k) {
               C.R_(i, j) += A.R_(i, k)*B.R_(k, j);
            }
         }
      }
 
      // Compute translation vector
      for (i = 0; i < D; ++i) {
         C.t_[i] = A.t_[i];
      }
      for (i = 0; i < D; ++i) {
         for (j = 0; j < D; ++j) {
            C.t_[i] += A.R_(i, j)*B.t_[j];
         }
      }
      C.normalize();
 
      return C;
   }
 
   /*
   * Matrix / IntVec<D> multiplication.
   */
   template <int D>
   IntVec<D> operator * (const SpaceSymmetry<D>& S, const IntVec<D>& V)
   {
      IntVec<D> U;
      int i, j;
      for (i = 0; i < D; ++i) {
         U[i] = 0;
         for (j = 0; j < D; ++j) {
            U[i] += S.R_(i,j)*V[j];
         }
      }
      return U;
   }
 
   /*
   * IntVec<D> / Matrix multiplication.
   */
   template <int D>
   IntVec<D> operator * (const IntVec<D>& V, const SpaceSymmetry<D>& S)
   {
      IntVec<D> U;
      int i, j;
      for (i = 0; i < D; ++i) {
         U[i] = 0;
         for (j = 0; j < D; ++j) {
            U[i] += V[j]*S.R_(j,i);
         }
      }
      return U;
   }
 
   /*
   * Output stream inserter for a SpaceSymmetry<D>.
   */ 
   template <int D>
   std::ostream& operator << (std::ostream& out, const SpaceSymmetry<D>& A)
   {
      int i, j;
      for (i = 0; i < D; ++i) {
         for (j = 0; j < D; ++j) {
            out << " " << Int(A.R_(i,j),2);
         }
         out << std::endl;
      }
      for (i = 0; i < D; ++i) {
         out << "  " << A.t_[i];
      }
      out << std::endl;
      return out;
   }
 
   /*
   * Input stream extractor for a SpaceSymmetry<D>.
   */ 
   template <int D>
   std::istream& operator >> (std::istream& in, SpaceSymmetry<D>& A)
   {
      int i, j;
      for (i = 0; i < D; ++i) {
         for (j = 0; j < D; ++j) {
            in >> A.R_(i,j);
         }
      }
      for (i = 0; i < D; ++i) {
         in >> A.t_[i];
      }
      return in;
   }
 
   // Define static members
   template<int D> 
   SpaceSymmetry<D> SpaceSymmetry<D>::identity_ = SpaceSymmetry<D>();
 
   template<int D> 
   bool SpaceSymmetry<D>::hasIdentity_ = false;
 
 
   #ifndef PRDC_SPACE_SYMMETRY_TPP
   // Suppress implicit instantiation
   extern template class SpaceSymmetry<1>;
   extern template class SpaceSymmetry<2>;
   extern template class SpaceSymmetry<3>;
   #endif 
 
}
}
#endif