SplitMonopole.H

Go to the documentation of this file.

 
#ifndef _SPLITMONOPOLE_H_
#define _SPLITMONOPOLE_H_
 
#include "Gas.H"
#include "../IO/ParmParser.H"
 
namespace criptic {
  namespace gas {
 
    class SplitMonopole : public Gas {
 
    public:
 
      SplitMonopole(const ParmParser& pp) {
 
    // Read data from ParmParser
    pp.get("gas.r0", r0);
    pp.get("gas.rho0", rho0);
    pp.get("gas.krho", krho);
    pp.get("gas.chi", chi);
    pp.get("gas.B0", B0);
    rFlat = 0.0;
    pp.query("gas.rflat", rFlat);
 
      }
 
      virtual ~SplitMonopole() { };
 
      virtual GasData gasData(const RealVec &x,
                  const Real t) const {
    GasData gd;
    Real r = x.mag();
    Real rfac = r > rFlat ? r / r0 : rFlat / r0;
    RealVec rhat = x / r;
 
    // Set gas properties
    gd.den = rho0 * pow(rfac, krho);
    gd.ionDen = gd.den * chi;
    gd.dx = r0;
    gd.v = zeroVec;
    gd.B = B0 / (rfac * rfac) * rhat;
    if (x[2] < 0) gd.B = -gd.B;
    gd.setAtomicComposition();
 
    // Set gradient properties; note that B = B0 (r/r0)^-2
    // \hat{r}, and with a bit of algebra one can show that
    // dB_i / dx_j = B / r^3 (r^2 delta_ij - 3 x_i x_j)
    if (r > rFlat)
      gd.denGrad = rho0 * krho * pow(rfac, krho) * rhat / r;
    else
      gd.denGrad = zeroVec;
    gd.ionDenGrad = gd.denGrad * chi;
    gd.vGrad = zeroTensor;
    if (r > rFlat)
      gd.BGrad = gd.B.mag() *
        (identityTensor/r - 3 * outer(x, x) / (r*r*r));
    else
      gd.BGrad = zeroTensor;
    if (x[2] < 0) gd.BGrad = -gd.BGrad;
 
    // Return
    return gd;
      }
      
      virtual void frame(const RealVec &x,
             const Real t,
             RealVec& v,
             TNBBasis& tnb) const {
    v = 0.0;
    tnb.eT = x.unit();
    if (x[2] < 0) tnb.eT = -tnb.eT;
    if (tnb.eT[0] == 0 && tnb.eT[1] == 0) {
      tnb.eN = xHat;
      tnb.eB = yHat;
    } else {
      tnb.eN = tnb.eT.cross(xHat).unit();
      tnb.eB = tnb.eT.cross(tnb.eN);
    }
      }
 
    private:
 
      Real r0;        
      Real rFlat;     
      Real rho0;      
      Real krho;      
      Real chi;       
      Real B0;        
    };
  }
}
#endif
// _SPLITMONOPOLE_H_