sphereCell.hh

#ifndef SPHERECELL_HH
#define SPHERECELL_HH
 
#include "sphereTopology.hh"
#include "geometry/cell2D.hh"
#include <cmath>
#include <random>
 
 
#define CONCEPTS_SPHERECELL_D 1
 
 
 
namespace concepts {
    
  //************************************************ SphereMapping ***
    
  struct SphereMapping {
    Real radius_;
    Real3d center_;
  };
    
  //************************************************ SphericalSurface3d ***
  
  class SphericalSurface3d : public Cell {
  public:
    SphericalSurface3d( SphericalSurface& cntr,
                        Real radius,
                        Real3d center );
      
    ~SphericalSurface3d() override;
    
    const Quad2d* child(uint i) const override;
    Quad2d* child(uint i) override;
    
    inline bool hasChildren() const { return chld_.size() > 0; }
 
    SphericalSurface& connector() const  override { return cntr_;}
    
    inline Real3d elemMap( const Real coord_local) const  override {
      throw concepts::MissingFeature("elementMap with one-dimensional coord cannot be called.");
    }
    
    inline Real3d elemMap( const Real2d& coord_local ) const  override {
      return elemMap(coord_local[0],coord_local[1]);
    }
    inline Real3d elemMap( const Real3d& coord_local ) const  override {
      return elemMap(coord_local[0],coord_local[1]);
    }
        
    inline Real3d elemMap( const Real& theta, const Real& phi ) const {
      return mapping_.center_ + mapping_.radius_*
        Real3d( std::cos(phi)*std::sin(theta),
                std::sin(phi)*std::sin(theta),
                std::cos(theta) ); 
    }    
    
    inline Real2d inverseElemMap( const Real3d& coord_global ) const {
      Real3d sphericalCoord = 1./mapping_.radius_*(coord_global-mapping_.center_);
      conceptsAssert( std::abs(sphericalCoord.l2()-1) < 1e3*std::numeric_limits<Real>::epsilon(),
                      Assertion() );
      Real phi = std::atan2(sphericalCoord[1],sphericalCoord[0]);
      if(phi < 0.) phi += 2*3.14159265358979323846264338328;
      Real theta = 0.;
      if( std::abs(sphericalCoord[2]) < 0.5 )
        theta =  std::acos(sphericalCoord[2]);     
      else {
        theta = std::asin( (std::abs(sphericalCoord[0]) > std::abs(sphericalCoord[1]) ? 
                               sphericalCoord[0]/std::cos(phi) : 
                               sphericalCoord[1]/std::sin(phi)) );
        if(sphericalCoord[2] < 0.)
          theta = 3.14159265358979323846264338328 - theta;
      }
      if(  std::abs(sphericalCoord[2]-1.) < 1e1*std::numeric_limits<concepts::Real>::epsilon()
        || std::abs(sphericalCoord[2]+1.) < 1e1*std::numeric_limits<concepts::Real>::epsilon()) {
        static std::default_random_engine engine_;
        static std::uniform_real_distribution<concepts::Real> distribution_(0.,2.*3.14159265358979323846264338328);
        phi = distribution_(engine_);
      }
      if( std::isnan( theta ) ) {
        if(std::abs(sphericalCoord[2]-1.)<1e1*std::numeric_limits<concepts::Real>::epsilon())
          theta = 1e1*std::numeric_limits<concepts::Real>::epsilon();
        else if(std::abs(sphericalCoord[2]+1.)<1e1*std::numeric_limits<concepts::Real>::epsilon())
          theta = 3.14159265358979323846264338328-1e1*std::numeric_limits<concepts::Real>::epsilon();
      }
      conceptsAssert(phi >= 0. && phi <= 2*3.14159265358979323846264338328,Assertion());      
      conceptsAssert(theta >= 0. && theta <= 3.14159265358979323846264338328,Assertion());
      conceptsAssert(!std::isnan(phi) && !std::isnan(theta),Assertion());      
      return Real2d( theta, phi/*std::acos(sphericalCoord[2]),
                     std::atan2(sphericalCoord[1],sphericalCoord[0])*/ );
    }
      
    inline Real jacobianDeterminant( const Real2d& coord_local ) const {
      return mapping_.radius_*mapping_.radius_*std::sin(coord_local[0]); }
      
    Real getRadius() const { return mapping_.radius_; }
    
    Real3d getCenter() const { return mapping_.center_; }    
    
    std::ostream& info(std::ostream& os) const override;
  private:
    SphericalSurface& cntr_;
    
    std::vector<Quad2d*> chld_;
    
    SphereMapping mapping_;
        
  };
  
  
  //************************************************************* Sphere3d ***
  
  class Sphere3d : public Cell {
  public:
    Sphere3d( Sphere& cntr,
              Real radius,
              Real3d center );
      
    ~Sphere3d() override;
    
    const Cell* child(uint i) const override;
    Quad2d* child(uint i) override;
    
    inline bool hasChildren() const { return false; }
 
    Sphere& connector() const  override { return cntr_;}
    
    inline Real3d elemMap( const Real coord_local ) const  override {
      throw concepts::MissingFeature("elementMap with one-dimensional coord cannot be called.");
    }
 
    inline Real3d elemMap( const Real2d& coord_local ) const  override {
      throw concepts::MissingFeature("elementMap with two-dimensional coord cannot be called.");
    }
 
    inline Real3d elemMap( const Real3d& coord_local ) const  override {
      return elemMap(coord_local[0],coord_local[1],coord_local[2]);
    }
    
    inline Real3d elemMap( const Real& radius, const Real& theta, const Real& phi ) const {
      return mapping_.center_ + mapping_.radius_*radius*
        Real3d( std::cos(phi)*std::sin(theta),
                std::sin(phi)*std::sin(theta),
                std::cos(theta) ); 
    }    
    
    inline Real3d inverseElemMap( const Real3d& coord_global ) const {
      Real3d sphericalCoord = 1./mapping_.radius_*(coord_global-mapping_.center_);
      const Real radius = sphericalCoord.l2();
      conceptsAssert(!std::isnan(radius),Assertion());
      conceptsAssert( std::abs(sphericalCoord.l2()-radius) < 1e2*std::numeric_limits<Real>::epsilon(),
                      Assertion() );
      sphericalCoord *= 1./radius;
      Real phi = std::atan2(sphericalCoord[1],sphericalCoord[0]);
      if(std::isnan(phi)) DEBUGL(1,"Phi is nan for spherical coordinates (" << sphericalCoord[0] << ","
                                  << sphericalCoord[1] << "," << sphericalCoord[2] << ")");
      if(phi < 0.) phi += 2*3.14159265358979323846264338328;
      Real theta = 0.;
      if( std::abs(sphericalCoord[2]) < 0.5 )
        theta =  std::acos(sphericalCoord[2]);     
      else {
        theta = std::asin( (std::abs(sphericalCoord[0]) > std::abs(sphericalCoord[1]) ? 
                               sphericalCoord[0]/std::cos(phi) : 
                               sphericalCoord[1]/std::sin(phi)) );
        if(sphericalCoord[2] < 0.)
          theta = 3.14159265358979323846264338328 - theta;
      }
      if(  std::abs(sphericalCoord[2]-1.) < 1e1*std::numeric_limits<concepts::Real>::epsilon()
        || std::abs(sphericalCoord[2]+1.) < 1e1*std::numeric_limits<concepts::Real>::epsilon()) {
        static std::default_random_engine engine_;
        static std::uniform_real_distribution<concepts::Real> distribution_(0.,2.*3.14159265358979323846264338328);
        phi = distribution_(engine_);
        conceptsAssert(phi <= 2.*3.14159265358979323846264338328,Assertion());
      }
      if( std::isnan( theta ) ) {
        if(std::abs(sphericalCoord[2]-1.)<1e1*std::numeric_limits<concepts::Real>::epsilon())
          theta = 1e1*std::numeric_limits<concepts::Real>::epsilon();
        else if(std::abs(sphericalCoord[2]+1.)<1e1*std::numeric_limits<concepts::Real>::epsilon())
          theta = 3.14159265358979323846264338328-1e1*std::numeric_limits<concepts::Real>::epsilon();
      }
      if ( !(phi >= 0. && phi <= 2*3.14159265358979323846264338328) ) {
        DEBUGL(1,"Phi value out of range: " << phi << ".");
        conceptsAssert(false,Assertion());
      }
      conceptsAssert(phi >= 0. && phi <= 2*3.14159265358979323846264338328,Assertion());      
      conceptsAssert(theta >= 0. && theta <= 2.*3.14159265358979323846264338328,Assertion());
      conceptsAssert(!std::isnan(phi) && !std::isnan(theta),Assertion());
      return Real3d( radius,
                     theta,
                     phi );
    }    
    
    inline Real jacobianDeterminant( const Real3d& coord_local ) const {
      return mapping_.radius_*mapping_.radius_*
             coord_local[0]*coord_local[0]*std::sin(coord_local[1]); }
      
    Real getRadius() const { return mapping_.radius_; }
    
    Real3d getCenter() const { return mapping_.center_; }    
      
    std::ostream& info(std::ostream& os) const override;
    
  private:
    Sphere& cntr_;
    
    SphereMapping mapping_;
    
  };
  
}
 
#endif // SPHERECELL_HH