Line data    Source code

       1             : #ifndef SET_BASE_H
       2             : #define SET_BASE_H
       3             : 
       4             : #include "AMReX.H"
       5             : #include "AMReX_REAL.H"
       6             : #include "AMReX_SPACE.H"
       7             : #include "AMReX_Vector.H"
       8             : #include "AMReX_BLassert.H"
       9             : #include "AMReX_Geometry.H"
      10             : #include "AMReX_IndexType.H"
      11             : #include <eigen3/Eigen/Core>
      12             : #include <eigen3/Eigen/Geometry>
      13             : #include <eigen3/Eigen/SVD>
      14             : #include <AMReX_REAL.H>
      15             : #include "Util/Util.H"
      16             : 
      17             : namespace Set
      18             : {
      19             : using Scalar   = amrex::Real;
      20             : using Vector   = Eigen::Matrix<amrex::Real,AMREX_SPACEDIM,1>;
      21             : using Vector3d = Eigen::Vector3d;
      22             : using Covector = Eigen::Matrix<amrex::Real,1,AMREX_SPACEDIM>;
      23             : using Matrix   = Eigen::Matrix<amrex::Real,AMREX_SPACEDIM,AMREX_SPACEDIM>;
      24             : using Matrix3d = Eigen::Matrix3d;
      25             : using iMatrix  = Eigen::Matrix<int,AMREX_SPACEDIM,AMREX_SPACEDIM>;
      26             : 
      27             : AMREX_FORCE_INLINE
      28             : Set::Matrix reduce(const Set::Matrix3d &in)
      29             : {
      30           0 :     return in.block<AMREX_SPACEDIM,AMREX_SPACEDIM>(0,0);
      31             : }
      32             : 
      33             : AMREX_FORCE_INLINE
      34             : Set::Matrix3d expand(const Set::Matrix &in)
      35             : {
      36           0 :     Set::Matrix3d ret = Set::Matrix3d::Zero();
      37           0 :     ret.block<AMREX_SPACEDIM,AMREX_SPACEDIM>(0,0) = in;
      38           0 :     return ret;
      39             : }
      40             : 
      41             : //using Quaternion = Eigen::Quaterniond;
      42             : 
      43             : class Quaternion : public Eigen::Quaterniond
      44             : {
      45             : public:
      46           8 :     Quaternion()
      47           8 :         : Eigen::Quaterniond() {}
      48         144 :     Quaternion(Set::Scalar w, Set::Scalar x, Set::Scalar y, Set::Scalar z)
      49         144 :         : Eigen::Quaterniond(w,x,y,z) {}
      50           0 :     Quaternion(const Eigen::Matrix3d &R)
      51           0 :         : Eigen::Quaterniond(R) {}
      52             :     Quaternion(Eigen::Quaterniond &q)
      53             :         : Eigen::Quaterniond(q.w(),q.x(),q.y(),q.z()) {}
      54             : 
      55             :     AMREX_FORCE_INLINE
      56             :     void operator = (const Eigen::Quaterniond &rhs)
      57             :     {
      58         128 :         w() = rhs.w(); x() = rhs.x(); y() = rhs.y(); z() = rhs.z();
      59         128 :     }
      60             : 
      61             :     AMREX_FORCE_INLINE
      62             :     void operator += (const Quaternion &rhs)
      63             :     {
      64           0 :         w() += rhs.w(); x() += rhs.x(); y() += rhs.y(); z() += rhs.z();
      65           0 :     }
      66             :     // AMREX_FORCE_INLINE
      67             :     // void operator * (const Set::Scalar alpha)
      68             :     // {
      69             :     //     w += alpha; x *= alpha; y *= alpha; z *= alpha;
      70             :     // }
      71             :     friend Quaternion operator * (const Set::Scalar alpha, const Quaternion b);
      72             :     friend Quaternion operator * (const Quaternion b,const Set::Scalar alpha);
      73             :     friend Quaternion operator + (const Quaternion a, const Quaternion b);
      74             :     friend Quaternion operator - (const Quaternion a, const Quaternion b);
      75             :     friend bool operator == (const Quaternion a, const Quaternion b);
      76             : };
      77             : AMREX_FORCE_INLINE
      78             : Quaternion operator * (const Set::Scalar alpha, const Quaternion b)
      79             : {
      80             :     Quaternion ret;
      81             :     ret.w() = b.w()*alpha; ret.x() = b.x()*alpha; ret.y() = b.y()*alpha; ret.z() = b.z()*alpha;
      82             :     return ret;
      83             : }    
      84             : AMREX_FORCE_INLINE
      85             : Quaternion operator * (const Quaternion b, const Set::Scalar alpha)
      86             : {
      87           8 :     Quaternion ret;
      88           8 :     ret.w() = b.w()*alpha; ret.x() = b.x()*alpha; ret.y() = b.y()*alpha; ret.z() = b.z()*alpha;
      89           8 :     return ret;
      90             : }    
      91             : AMREX_FORCE_INLINE
      92             : Quaternion operator + (const Quaternion a, const Quaternion b)
      93             : {
      94           0 :     Quaternion ret;
      95           0 :     ret.w() = a.w()+b.w(); ret.x() = a.x()+b.x(); ret.y() = a.y()+b.y(); ret.z() = a.z()+b.z();
      96           0 :     return ret;
      97             : }    
      98             : 
      99             : AMREX_FORCE_INLINE
     100             : Quaternion operator - (const Quaternion a, const Quaternion b)
     101             : {
     102           0 :     Quaternion ret;
     103           0 :     ret.w() = a.w()-b.w(); ret.x() = a.x()-b.x(); ret.y() = a.y()-b.y(); ret.z() = a.z()-b.z();
     104           0 :     return ret;
     105             : }    
     106             : 
     107             : AMREX_FORCE_INLINE
     108             : bool operator == (const Quaternion a, const Quaternion b)
     109             : {
     110          20 :     if (fabs((a.w() - b.w())/(a.w()+a.w())) > 1E-8) return false;
     111          20 :     if (fabs((a.x() - b.x())/(a.x()+a.x())) > 1E-8) return false;
     112          20 :     if (fabs((a.y() - b.y())/(a.y()+a.y())) > 1E-8) return false;
     113          20 :     if (fabs((a.z() - b.z())/(a.z()+a.z())) > 1E-8) return false;
     114          20 :     return true;
     115             : }
     116             : 
     117             : 
     118             : static Scalar Garbage = NAN;
     119             : 
     120             : AMREX_FORCE_INLINE
     121             : Vector Position(const int &i, const int &j, const int &k, const amrex::Geometry &geom, const amrex::IndexType &ixType)
     122             : {
     123             :     (void)k;
     124      204836 :     if (ixType == amrex::IndexType::TheNodeType())
     125             :     {
     126      146596 :         return Vector(AMREX_D_DECL(
     127             :                 geom.ProbLo()[0] + ((amrex::Real)(i)) * geom.CellSize()[0],
     128             :                 geom.ProbLo()[1] + ((amrex::Real)(j)) * geom.CellSize()[1],
     129       73298 :                 geom.ProbLo()[2] + ((amrex::Real)(k)) * geom.CellSize()[2]));
     130             :     }
     131       58240 :     else if (ixType == amrex::IndexType::TheCellType())
     132             :     {
     133       58240 :         return Vector(AMREX_D_DECL(
     134             :                 geom.ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom.CellSize()[0],
     135             :                 geom.ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom.CellSize()[1],
     136       29120 :                 geom.ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom.CellSize()[2]));
     137             :     }
     138           0 :     Util::Abort(INFO,"Unsupported index type");
     139           0 :     return Set::Vector::Zero(); 
     140             : }
     141             : 
     142             : AMREX_FORCE_INLINE
     143             : Vector Normal(  AMREX_D_DECL(bool xmin,bool ymin, bool zmin),
     144             :                 AMREX_D_DECL(bool xmax,bool ymax, bool zmax))
     145             : {
     146           0 :     Set::Vector N =  Set::Vector::Zero();
     147           0 :     if (xmin)      N(0) = -1.0;
     148           0 :     else if (xmax) N(0) =  1.0;
     149             : #if AMREX_SPACEDIM>1
     150           0 :     if (ymin)      N(1) = -1.0;
     151           0 :     else if (ymax) N(1) =  1.0;
     152             : #endif
     153             : #if AMREX_SPACEDIM>2
     154           0 :     if (zmin)      N(2) = -1.0;
     155           0 :     else if (zmax) N(2) =  1.0;
     156             : #endif 
     157           0 :     return N;
     158             : }
     159             : 
     160             : AMREX_FORCE_INLINE
     161             : Vector Size(const amrex::Geometry &geom)
     162             : {
     163           0 :     Set::Vector size;
     164           0 :     size(0) = geom.ProbHi()[0] - geom.ProbLo()[0];
     165             : #if AMREX_SPACEDIM>1
     166           0 :     size(1) = geom.ProbHi()[1] - geom.ProbLo()[1];
     167             : #endif
     168             : #if AMREX_SPACEDIM>2
     169           0 :     size(2) = geom.ProbHi()[2] - geom.ProbLo()[2];
     170             : #endif
     171           0 :     return size;
     172             : }
     173             : 
     174             : 
     175             : AMREX_FORCE_INLINE
     176             : Vector Volume(const int &i, const int &j, const int &k, const amrex::Geometry &geom, const amrex::IndexType &ixType)
     177             : {
     178             :     (void)k;
     179             :     if (ixType == amrex::IndexType::TheNodeType())
     180             :     {
     181             :         return Vector(AMREX_D_DECL(
     182             :                 geom.ProbLo()[0] + ((amrex::Real)(i)) * geom.CellSize()[0],
     183             :                 geom.ProbLo()[1] + ((amrex::Real)(j)) * geom.CellSize()[1],
     184             :                 geom.ProbLo()[2] + ((amrex::Real)(k)) * geom.CellSize()[2]));
     185             :     }
     186             :     else if (ixType == amrex::IndexType::TheCellType())
     187             :     {
     188             :         return Vector(AMREX_D_DECL(
     189             :                 geom.ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom.CellSize()[0],
     190             :                 geom.ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom.CellSize()[1],
     191             :                 geom.ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom.CellSize()[2]));
     192             :     }
     193             :     Util::Abort(INFO,"Unsupported index type");
     194             :     return Set::Vector::Zero(); 
     195             : }
     196             : 
     197             : enum Sym {None,Major,Minor,MajorMinor,Diagonal,Full,Isotropic};
     198             : template<int dim,int sym> class Matrix4{};
     199             : template<int dim, int sym>
     200             : std::ostream&
     201           7 : operator<< (std::ostream& os, const Matrix4<dim,sym>& b)
     202             : {
     203           7 :     Matrix4<dim,sym> bcopy = b;
     204           7 :         bcopy.Print(os);
     205           7 :     return os;
     206             : }
     207             : }
     208             : #endif