/* title
   在二維網格區域內畫等值線
 */

/* comment *****************************************************
   
  讀入三角或四邊形網格的資料與點的數值，程式可畫出等值線



   由 fem.dat 檔案讀入資料後，以 gnuplot 資料型式輸出畫等值線

   檔案格式：

     1 ： 輸出檔的名稱

     2 ： 等值線的最小值  等值線的最大值  等值線的等份數  

     3 ： 每個格式的切割等份數 ( 數字越大，等值線越平滑 )

     4 ： 格子點的多項式插分階數 ( 1 或 2 )

     5 ： 輸入(三角形或四邊形)格子的所有點的座標與其對應數值

   說明：
     (1) 資料檔的檔名為 fem.dat
     (2) 每列的第一個字元若為井號，則此列為註解列
     (3) 不管是三角形或四邊形，格子點的順序依次為
         頂點，邊的中點，內部點
     (4) 格子與格子之間的資料以空行分開
     (5) 各格子的點數為：

               階數 |  三角形    四邊形
              ==========================
                 1  |     3         4
                 2  |     6         9

     (6) gnuplot 內使用以下指令看等值線 
            plot "輸出檔檔名" with line 


 ******************************************************
   以下這裡所用的 fem.dat 資料檔格式範例：
 ******************************************************

        # 輸出檔名稱
        aa.dat
        
        # 等值線的最小值，最大值，與等份數
        0 100 40 
        
        # 每個格子的切割數 (數字越大，等值線越平滑)
        40
        
        # 格子的插分多項數的階數
        2
        
        # 三角形
        0  0  30
        10 0  20
        0 10  15
        5  0  20
        5  5  17
        0  5   5
        
        # 三角形
        10 0  20
        10 10 30
        0  10 15
        10 5  25
        5 10  10
        5 5   17
        
        # 四邊形
        10 0  20
        20 0  30
        20 10 30
        10 10 30
        15 0  30 
        20 5  35
        15 10 50
        10 5  25
        15 5  40
                
******************************************************/

#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include <string>
#include <iterator>
#include <algorithm>
#include <map>

using namespace std ;

// vector point class
class  Point {
  private :
    double  x , y ;

  public :

    Point( double a = 0. , double b = 0. ) : x(a) , y(b) {}

    double  getx() const { return x ; }
    double  gety() const { return y ; }

    void    setx( double a ) { x = a ; }
    void    sety( double a ) { y = a ; }

    void    set_point( double a , double b ) { x = a ; y = b ; }

    
    friend  istream& operator>> ( istream& in , Point& pt ) {
        return  in >> pt.x >> pt.y ;
    }

    friend  ostream& operator<< ( ostream& out , const Point& pt ) {
        return  out << pt.x << " " << pt.y ;
    }
    
};


bool  operator== ( const Point& a , const Point& b ) {
    return ( a.getx() == b.getx() && a.gety() == b.gety() ) ;
}


bool  operator< ( const Point& a , const Point& b ) {
    if ( a.getx() != b.getx() ) 
        return  a.getx() < b.getx() ;
    else 
        return  a.gety() < b.gety() ;
}


Point  operator*( double c , const Point& pt ) {
    return  Point(c*pt.getx(),c*pt.gety()) ;
}

Point  operator*( const Point& pt , double c ) {
    return  Point(c*pt.getx(),c*pt.gety()) ;
}

Point  operator/( const Point& pt , double c ) {
    return  Point(pt.getx()/c,pt.gety()/c) ;
}

Point  operator+( const Point& a , const Point& b ) {
    return  Point(a.getx()+b.getx(),a.gety()+b.gety()) ;
}

Point  operator-( const Point& a , const Point& b ) {
    return  Point(a.getx()-b.getx(),a.gety()-b.gety()) ;
}

    


// 抽象基礎類別 ： shape function 
struct Shape_Function {
    virtual  double  value( const Point& ) const = 0 ;
    virtual ~Shape_Function() = 0 ;
};
Shape_Function::~Shape_Function() {}


// 矩形格子形狀函數
struct Rectangular_Shape_Function : public Shape_Function {
    virtual  double  value( const Point& ) const = 0 ;
    virtual ~Rectangular_Shape_Function() = 0 ;
};
Rectangular_Shape_Function::~Rectangular_Shape_Function() {}


// 一階矩形格子形狀函數
struct First_Order_RSF : public Rectangular_Shape_Function {
    virtual  double  value( const Point& ) const = 0 ;
    virtual ~First_Order_RSF() = 0 ;
};
First_Order_RSF::~First_Order_RSF() {}


// first function for first order shape function
struct RSF11 : public First_Order_RSF {
    double  value( const Point& pt ) const {
        return  0.25*(1.-pt.getx())*(1.-pt.gety()) ;
    }
};
    
// second function for first order shape function
struct RSF12 : public First_Order_RSF {
    double  value( const Point& pt ) const {
        return  0.25*(1.+pt.getx())*(1.-pt.gety()) ;
    }
};

// third function for first order shape function
struct RSF13 : public First_Order_RSF {
    double  value( const Point& pt ) const {
        return  0.25*(1.+pt.getx())*(1.+pt.gety()) ;
    }
};

// fourth function for first order shape function
struct RSF14 : public First_Order_RSF {
    double  value( const Point& pt ) const {
        return  0.25*(1.-pt.getx())*(1.+pt.gety()) ;
    }
};


// 二階矩形格子形狀函數
struct Second_Order_RSF : public Rectangular_Shape_Function {
    virtual  double  value ( const Point& ) const = 0 ;
    virtual ~Second_Order_RSF() = 0 ;
};
Second_Order_RSF::~Second_Order_RSF() {}

// first function for second order shape function
struct RSF21 : public Second_Order_RSF {
    double  value( const Point& pt ) const {
        return  0.25*(1.-pt.getx())*pt.getx()*(1.-pt.gety())*pt.gety() ;
    }
};
    
// second function for second order shape function
struct RSF22 : public Second_Order_RSF {
    double  value( const Point& pt ) const {
        return  -0.25*(1.+pt.getx())*pt.getx()*(1.-pt.gety())*pt.gety() ;
    }
};

// third function for second order shape function
struct RSF23 : public Second_Order_RSF {
    double  value( const Point& pt ) const {
        return  0.25*(1.+pt.getx())*pt.getx()*(1.+pt.gety())*pt.gety() ;
    }
};

// fourth function for second order shape function
struct RSF24 : public Second_Order_RSF {
    double  value( const Point& pt ) const {
        return  -0.25*(1.-pt.getx())*pt.getx()*(1.+pt.gety())*pt.gety() ;
    }
};

// fifth function for second order shape function
struct RSF25 : public Second_Order_RSF {
    double  value( const Point& pt ) const {
        return  -0.5*(1.-pt.getx())*(1.+pt.getx())*(1.-pt.gety())*pt.gety() ;
    }
};
    
// sixth function for second order shape function
struct RSF26 : public Second_Order_RSF {
    double  value( const Point& pt ) const {
        return  0.5*(1.+pt.getx())*pt.getx()*(1.-pt.gety())*(1.+pt.gety()) ;
    }
};

// seventh function for second order shape function
struct RSF27 : public Second_Order_RSF {
    double  value( const Point& pt ) const {
        return  0.5*(1.-pt.getx())*(1.+pt.getx())*(1.+pt.gety())*pt.gety() ;
    }
};

// eighth function for second order shape function
struct RSF28 : public Second_Order_RSF {
    double  value( const Point& pt ) const {
        return  -0.5*(1.-pt.getx())*pt.getx()*(1.+pt.gety())*(1.-pt.gety()) ;
    }
};

// nineth function for second order shape function
struct RSF29 : public Second_Order_RSF {
    double  value( const Point& pt ) const {
        return  (1.-pt.getx())*(1.+pt.getx())*(1.+pt.gety())*(1.-pt.gety()) ;
    }
};





// 三角形格子形狀函數
struct Triangular_Shape_Function : public Shape_Function {
    virtual  double  value( const Point& ) const = 0 ;
    virtual ~Triangular_Shape_Function() = 0 ;
};
Triangular_Shape_Function::~Triangular_Shape_Function() {}


// 一階三角形格子形狀函數
struct First_Order_TSF : public Triangular_Shape_Function {
    virtual  double  value( const Point& ) const = 0 ;
    virtual ~First_Order_TSF() = 0 ;
};
First_Order_TSF::~First_Order_TSF() {}


// first function for first order shape function
struct TSF11 : public First_Order_TSF {
    double  value( const Point& pt ) const {
        return  (1.-pt.getx()-pt.gety()) ;
    }
};
    
// second function for first order shape function
struct TSF12 : public First_Order_TSF {
    double  value( const Point& pt ) const {
        return  pt.getx() ;
    }
};

// third function for first order shape function
struct TSF13 : public First_Order_TSF {
    double  value( const Point& pt ) const {
        return  pt.gety() ;
    }
};



// 二階三角形格子形狀函數
struct Second_Order_TSF : public Triangular_Shape_Function {
    virtual  double  value( const Point& ) const = 0 ;
    virtual ~Second_Order_TSF() = 0 ;
};
Second_Order_TSF::~Second_Order_TSF() {}


// first function for second order shape function
struct TSF21 : public Second_Order_TSF {
    double  value( const Point& pt ) const {
        return  2.*(0.5-pt.getx()-pt.gety())*(1-pt.getx()-pt.gety()) ;
    }
};

// second function for second order shape function
struct TSF22 : public Second_Order_TSF {
    double  value( const Point& pt ) const {
        return  -2*pt.getx()*(0.5-pt.getx()) ;
    }
};

// third function for second order shape function
struct TSF23 : public Second_Order_TSF {
    double  value( const Point& pt ) const {
        return  -2*pt.gety()*(0.5-pt.gety()) ;
    }
};

// fourth function for second order shape function
struct TSF24 : public Second_Order_TSF {
    double  value( const Point& pt ) const {
        return  4*pt.getx()*(1-pt.getx()-pt.gety()) ;
    }
};

// fifth function for second order shape function
struct TSF25 : public Second_Order_TSF {
    double  value( const Point& pt ) const {
        return  4*pt.getx()*pt.gety() ;
    }
};

// sixth function for second order shape function
struct TSF26 : public Second_Order_TSF {
    double  value( const Point& pt ) const {
        return  4*pt.gety()*(1.-pt.getx()-pt.gety()) ;
    }
};
    


// FEM 形狀函數基礎類別
struct FEM_Shape_Function {
    int deg ;
    FEM_Shape_Function( int d ) : deg(d) {}
    virtual ~FEM_Shape_Function() = 0 ;
};
FEM_Shape_Function::~FEM_Shape_Function() {}


// FEM 矩形格子形狀函數類別
struct FEM_Rectangular_Shape_Function : public FEM_Shape_Function {

    FEM_Rectangular_Shape_Function( int d ) : FEM_Shape_Function(d) {}

    vector<Rectangular_Shape_Function*>  get_shape_fn() const {

        vector<Rectangular_Shape_Function*>  sf ;

        if ( deg == 1 ) {
            sf.push_back( new RSF11 ) ;
            sf.push_back( new RSF12 ) ;
            sf.push_back( new RSF13 ) ;
            sf.push_back( new RSF14 ) ;
        } else if ( deg == 2 ) {
            sf.push_back( new RSF21 ) ;
            sf.push_back( new RSF22 ) ;
            sf.push_back( new RSF23 ) ;
            sf.push_back( new RSF24 ) ;
            sf.push_back( new RSF25 ) ;
            sf.push_back( new RSF26 ) ;
            sf.push_back( new RSF27 ) ;
            sf.push_back( new RSF28 ) ;
            sf.push_back( new RSF29 ) ;
        }
        return  sf ;

    }
    
};
    

FEM_Rectangular_Shape_Function  R_First_Order(1) ;
FEM_Rectangular_Shape_Function  R_Second_Order(2) ;


// FEM 三角形格子形狀函數類別
struct FEM_Triangular_Shape_Function : public FEM_Shape_Function {

    FEM_Triangular_Shape_Function( int d ) : FEM_Shape_Function(d) {}

    vector<Triangular_Shape_Function*>  get_shape_fn() const {

        vector<Triangular_Shape_Function*>  sf ;

        if ( deg == 1 ) {
            sf.push_back( new TSF11 ) ;
            sf.push_back( new TSF12 ) ;
            sf.push_back( new TSF13 ) ;
        } else if ( deg == 2 ) {
            sf.push_back( new TSF21 ) ;
            sf.push_back( new TSF22 ) ;
            sf.push_back( new TSF23 ) ;
            sf.push_back( new TSF24 ) ;
            sf.push_back( new TSF25 ) ;
            sf.push_back( new TSF26 ) ;
        }
        
        return  sf ;

    }
    
};


FEM_Triangular_Shape_Function  T_First_Order(1) ;
FEM_Triangular_Shape_Function  T_Second_Order(2) ;



class  Cell {
  protected :
    int             no_pts ;        // no of points in cell
    int             deg ;           // interpolant polynomial of degree
    
    vector<Point>   pts ;           // points on cell
    vector<double>  vals ;          // value at each point


  public :

    Cell( const vector<Point>& p , const vector<double>& v ) : pts(p) , vals(v) {
        no_pts = pts.size() ;           // 儲存格子總點數
    }


    // 回傳 no of edge
    virtual  int  no_edge() const = 0 ;

    // 回傳格子的總點數
    virtual  int no_cell_pts() const = 0 ;

    // 若主元素 (master element) 的每邊切割成 m 個等份，則
    // 函式回傳在下標區域 (s,t) 頂點的實際位置
    virtual  vector<Point>  get_grid_points( int s , int t , int m ) const = 0 ;

    // 回傳在主座標系統下的 p 點數值插分值
    virtual double  get_val( const Point& p ) const = 0 ;

    // 回傳 每邊切割成 m 個等份的 (s,t) 區域的頂點的插分數值
    virtual vector<double>  get_grid_val( int s , int t , int m ) const = 0 ;

    // 輸入主座標系統的格子位置與其對應數值，
    // 函式回傳一等值 val 直線的兩個端點，
    // 此直線連接格子的第 s 邊到第 t 邊
    virtual vector<Point>  division_isoline( const vector<Point>& p , 
                                             const vector<double>& v ,
                                             int s , int t , double val ) const = 0 ;
    
    // 找出間格值陣列 isovalues 在 [a,b] 之間的元素陣列
    virtual vector<double>  find_subvalues( double a , double b , 
                                            const vector<double>& isovalues ) const = 0 ;

    // 回傳本格子的四個邊界點實際位置
    virtual vector<Point>  cell_border_pts() const = 0 ;

    // 將數值 [a,b] 區間設定 n 個等份， 輸出等數值限到檔案 fn 
    // Cell 的每邊切割成 m 個小等份，
    virtual void  plot_isogram( ostream& outfile , double a , double b , int n , int m ,
                                bool  plot_cell_boundary ) const = 0 ;
    

    virtual ~Cell() = 0 ;
        
};


Cell::~Cell() {}
 


// 方形格子插分
//
class  Rectangular_Cell : public Cell {
  private :

    // shape function 靜態資料成員
    static vector<Rectangular_Shape_Function*>  shape_fn ;

    enum  { No_Edge = 4 };

    
    // 根據 d 值設定靜態資料成員 shape_fn
    static  void  setup_shape_function( int d ) {
        if ( d == 1 ) 
            shape_fn = R_First_Order.get_shape_fn() ;
        else if ( d == 2 )
            shape_fn = R_Second_Order.get_shape_fn() ;
    }


    // 若主元素 (master element) 的每邊切割成 m 個等份，則
    // 函式回傳在下標區域 (s,t) 四個頂點的實際位置
    // s , t 為介於 [0,m-1] 的整數
    vector<Point>  get_grid_points( int s , int t , int m ) const {

        int      i , j ;
        double  xa , ya ;

        // p     : 主座標系統位置 (point in master coord)
        // s , t : 介於 [0,m-1]
        Point   p[No_Edge] = { Point(-1+s*2./m,-1+t*2./m) ,         
                               Point(-1+(s+1)*2./m,-1+t*2./m) , 
                               Point(-1+(s+1)*2./m,-1+(t+1)*2./m) , 
                               Point(-1+s*2./m,-1+(t+1)*2./m) };
        
        // q : p 所對應的實際座標系統位置 (point in phyical coord)
        vector<Point>  q(No_Edge) ;

        // 根據 p 點在主座標系統計算其在實際座標系統位置
        for ( i = 0 ; i < No_Edge ; ++i ) {
            xa = ya = 0. ;
            for ( j = 0 ; j < no_pts ; ++j ) {
                xa += pts[j].getx() * shape_fn[j]->value(p[i]) ;
                ya += pts[j].gety() * shape_fn[j]->value(p[i]) ;
            }
            q[i].set_point(xa,ya) ;
        }
        return  q ;
        
    }
    
    // 回傳在主座標系統下的 p 點數值插分值
    double  get_val( const Point& p ) const {
        double  v = 0 ;
        for ( int i = 0 ; i < no_pts ; ++i ) {
            v += vals[i] * shape_fn[i]->value(p) ;
        }
        return v ;
    }

    // 回傳 每邊切割成 m 個等份的 (s,t) 區域的頂點的插分數值
    vector<double>  get_grid_val( int s , int t , int m ) const {

        Point   p[No_Edge] = { Point(-1+s*2./m,-1+t*2./m) ,         
                               Point(-1+(s+1)*2./m,-1+t*2./m) , 
                               Point(-1+(s+1)*2./m,-1+(t+1)*2./m) , 
                               Point(-1+s*2./m,-1+(t+1)*2./m) } ;

        vector<double>  v(No_Edge) ;
        for ( int i = 0 ; i < No_Edge ; ++i ) v[i] = get_val(p[i]) ;
        return  v ;

    }

    
    // 輸入主座標系統的格子位置與其對應數值，
    // 函式回傳一等值 val 直線的兩個端點，
    // 此直線連接格子的第 s 邊到第 t 邊
    // s , t 介於 [0,3] 之間
    vector<Point>  division_isoline( 
        const vector<Point>& p , const vector<double>& v ,
        int s , int t , double val ) const {

        double  d ;
        Point   p1 , p2 , p3 , p4 ;
        double  v1 , v2 , v3 , v4 ;
        
        vector<Point>  line_pts ;
        
        p1 = p[s] ;
        v1 = v[s] ;

        p2 = p[(s+1)%No_Edge] ;
        v2 = v[(s+1)%No_Edge] ;

        // 如果 s 邊上兩端點的數值相等，且與設定數值相等
        if ( v1 == v2 && v1 == val ) {
            line_pts.push_back(p1) ;
            line_pts.push_back(p2) ;
            return  line_pts ;
        } else if ( v1 == v2 || v3 == v4 ) {
            return  line_pts ;
        }

        // 起點位置
        if ( v1 <= val && val < v2 ) {
            d = (val-v1)/(v2-v1) ;
            line_pts.push_back( p1 + d * (p2-p1) ) ;
        } else if ( v2 <= val && val < v1 ) {
            d = (val-v2)/(v1-v2) ;
            line_pts.push_back( p2 + d * (p1-p2) ) ;
        } else {
            return  line_pts ;
        }


        p3 = p[t] ;
        v3 = v[t] ;

        p4 = p[(t+1)%No_Edge] ;
        v4 = v[(t+1)%No_Edge] ;

        // 終點位置
        if ( v3 <= val && val < v4 ) {
            d = (val-v3)/(v4-v3) ;
            line_pts.push_back( p3 + d * (p4-p3) ) ;
        } else if ( v4 <= val && val < v3 ) {
            d = (val-v4)/(v3-v4) ;
            line_pts.push_back( p4 + d * (p3-p4) ) ;
        }
         
        
        return  line_pts ;

    }

    // 找出間格值陣列 isovalues 在 [a,b] 之間的元素陣列
    vector<double>  find_subvalues( double a , double b , 
                                    const vector<double>& isovalues ) const {
        vector<double>  foo ;

        double  max0 = ( a > b ? a : b ) ;
        double  min0 = ( a > b ? b : a ) ;
        
        for ( int i = 0 ; i < isovalues.size() ; ++i ) {
            if ( min0 <= isovalues[i] && isovalues[i] <= max0 ) 
                foo.push_back(isovalues[i]) ;
        }

        return  foo ;

    }
    
    
  public :

    // 建構函式
    Rectangular_Cell( const vector<Point>& p , const vector<double>& v ) 
        : Cell(p,v) {

        // 找出 deg
        deg = 1 ;
        while ( (deg+1)*(deg+1) < no_pts ) ++deg ;

        // 根據 deg 設定使用的 shape function
        if ( shape_fn.size() == 0 ) setup_shape_function(deg) ;

    }

    // 回傳 no of edge
    int  no_edge() const { return  static_cast<int>(No_Edge) ; }
    
    // 回傳格子的總點數
    int  no_cell_pts() const { return (deg+1)*(deg+1) ; }

    // 功能與上一樣，但設為 靜態函式
    static  int  no_cell_pts(int d) { return (d+1)*(d+1) ; }

    // 回傳本格子的四個邊界點實際位置
    vector<Point>  cell_border_pts() const {
        return  vector<Point>(pts.begin(),pts.begin()+No_Edge) ;
    }
    
        
    // 將數值 [a,b] 區間設定 n 個等份， 輸出等數值限到檔案 fn 
    // Cell 的每邊切割成 m 個小等份，
    void  plot_isogram( ostream& outfile , double a , double b , int n , int m ,
                        bool  plot_cell_boundary ) const {

        int             i , j , k , s , t ;
        double          dx = 2./m ;
        vector<Point>   division_pts , border_pts ;
        vector<double>  iso_values(n+1) ;
        vector<double>  division_vals ;

        vector<Point>   isoline ;
        vector<double>  sub_isovals ;


        // 計算所要畫出的所有等值點
        double          dv = (b-a)/n ;
        for ( i = 0 ; i <= n ; ++i ) iso_values[i] = a + i*dv ;


        // 畫格子的邊界
        if ( plot_cell_boundary ) {
            border_pts = cell_border_pts() ;
            copy(border_pts.begin(),border_pts.end(),ostream_iterator<Point>(outfile,"\n")) ;
            outfile << border_pts[0] << "\n\n" ;
        }
        

        // 將格子每邊各切割成 m 等份，然後畫出每一個小區間範圍的等值線
        //
        for ( i = 0 ; i < m ; ++i ) {

            for ( j = 0 ; j < m ; ++j ) {

                // 找出圍住 (i,j) 小區間四個頂點的實際位置與其對應數值
                division_pts  = get_grid_points(i,j,m) ;
                division_vals = get_grid_val(i,j,m) ;

                // 找出每個小區間的 s 邊與剩餘的 t 邊的等值線的兩個端點
                for ( s = 0 ; s < No_Edge ; ++s ) {

                    // 找出在小區間的 s 邊上的實際可用的等值陣列
                    sub_isovals = find_subvalues(division_vals[s],
                                                 division_vals[(s+1)%4],
                                                 iso_values) ;

                    // 找出 s 到 t 邊的等值線
                    for ( t = (s+1)%No_Edge ; t < No_Edge ; ++t ) {

                        for ( k = 0 ; k < sub_isovals.size() ; ++k ) {
                            isoline = division_isoline(division_pts,division_vals,
                                                       s,t,sub_isovals[k]) ;
                            if ( isoline.size() == 2 ) {
                                outfile << isoline[0] << "\n" << isoline[1] << "\n\n" ;
                            }
                        }

                        if ( s == No_Edge-1 ) break ;
                        
                    }
                }

            }
        }
        
    }

    ~Rectangular_Cell() {}
    
};


vector<Rectangular_Shape_Function*>  Rectangular_Cell::shape_fn ;



// 三角形格子插分
//
class  Triangular_Cell : public Cell {
  private :

    // shape function 靜態資料成員
    static vector<Triangular_Shape_Function*>  shape_fn ;

    enum { No_Edge = 3 } ;  


    // 根據 d 值設定靜態資料成員 shape_fn
    static  void  setup_shape_function( int d ) {
        if ( d == 1 ) 
            shape_fn = T_First_Order.get_shape_fn() ;
        else if ( d == 2 ) 
            shape_fn = T_Second_Order.get_shape_fn() ;
    }


    // 若主元素 (master element) 的每邊切割成 m 個等份，則
    // 函式回傳在下標區域 (s,t) 三個頂點的實際位置
    // s , t : [0,2(m-1)] ,   s+t in [0,2(m-1)]
    vector<Point>  get_grid_points( int s , int t , int m ) const {

        int      i , j ;
        double  xa , ya ;
        double  dh = 1./m ;

        // p     : 主座標系統位置 (point in master coord)
        // s , t : 介於 [0,2(m-1)]
        Point   p[No_Edge] ;

        xa = ( s%2 ? 0.5*(s+1)/m : 0.5*s/m ) ;
        ya = ( t%2 ? 0.5*(t+1)/m : 0.5*t/m ) ;

        p[0] = Point(xa,ya) ;
        
        if ( s % 2 == 0 ) {
            p[1] = Point(xa+dh,ya) ;
            p[2] = Point(xa,ya+dh) ;
        } else {
            p[1] = Point(xa-dh,ya) ;
            p[2] = Point(xa,ya-dh) ;
        }
        
        
        
        
        // q : p 所對應的實際座標系統位置 (point in phyical coord)
        vector<Point>  q(No_Edge) ;

        // 根據 p 點在主座標系統計算其在實際座標系統位置
        for ( i = 0 ; i < No_Edge ; ++i ) {
            xa = ya = 0. ;
            for ( j = 0 ; j < no_pts ; ++j ) {
                xa += pts[j].getx() * shape_fn[j]->value(p[i]) ;
                ya += pts[j].gety() * shape_fn[j]->value(p[i]) ;
            }
            q[i].set_point(xa,ya) ;
        }
        return  q ;
        
    }
    
    // 回傳在主座標系統下的 p 點數值插分值
    double  get_val( const Point& p ) const {
        double  v = 0 ;
        for ( int i = 0 ; i < no_pts ; ++i ) {
            v += vals[i] * shape_fn[i]->value(p) ;
        }
        return v ;
    }

    // 回傳 每邊切割成 m 個等份的 (s,t) 區域的頂點的插分數值
    vector<double>  get_grid_val( int s , int t , int m ) const {

        // p     : 主座標系統位置 (point in master coord)
        // s , t : 介於 [0,2(m-1)]
        Point   p[No_Edge] ;

        double  xa = ( s%2 ? 0.5*(s+1)/m : 0.5*s/m ) ;
        double  ya = ( t%2 ? 0.5*(t+1)/m : 0.5*t/m ) ;
        double  dh = 1./m ;

        p[0] = Point(xa,ya) ;
        
        if ( s % 2 == 0 ) {
            p[1] = Point(xa+dh,ya) ;
            p[2] = Point(xa,ya+dh) ;
        } else {
            p[1] = Point(xa-dh,ya) ;
            p[2] = Point(xa,ya-dh) ;
        }
        

        vector<double>  v(No_Edge) ;
        for ( int i = 0 ; i < No_Edge ; ++i ) v[i] = get_val(p[i]) ;
        return  v ;

    }

    
    // 輸入主座標系統的格子位置與其對應數值，
    // 函式回傳一等值 val 直線的兩個端點，
    // 此直線連接格子的第 s 邊到第 t 邊
    // s , t 介於 [0,2] 之間
    vector<Point>  division_isoline( 
        const vector<Point>& p , const vector<double>& v ,
        int s , int t , double val ) const {

        double  d ;
        Point   p1 , p2 , p3 , p4 ;
        double  v1 , v2 , v3 , v4 ;
        
        vector<Point>  line_pts ;
        
        p1 = p[s] ;
        v1 = v[s] ;

        p2 = p[(s+1)%No_Edge] ;
        v2 = v[(s+1)%No_Edge] ;

        // 如果 s 邊上兩端點的數值相等，且與設定數值相等
        if ( v1 == v2 && v1 == val ) {
            line_pts.push_back(p1) ;
            line_pts.push_back(p2) ;
            return  line_pts ;
        } else if ( v1 == v2 || v3 == v4 ) {
            return  line_pts ;
        }

        // 起點位置
        if ( v1 <= val && val < v2 ) {
            d = (val-v1)/(v2-v1) ;
            line_pts.push_back( p1 + d * (p2-p1) ) ;
        } else if ( v2 <= val && val < v1 ) {
            d = (val-v2)/(v1-v2) ;
            line_pts.push_back( p2 + d * (p1-p2) ) ;
        } else {
            return  line_pts ;
        }


        p3 = p[t] ;
        v3 = v[t] ;

        p4 = p[(t+1)%No_Edge] ;
        v4 = v[(t+1)%No_Edge] ;

        // 終點位置
        if ( v3 <= val && val < v4 ) {
            d = (val-v3)/(v4-v3) ;
            line_pts.push_back( p3 + d * (p4-p3) ) ;
        } else if ( v4 <= val && val < v3 ) {
            d = (val-v4)/(v3-v4) ;
            line_pts.push_back( p4 + d * (p3-p4) ) ;
        }
         
        
        return  line_pts ;

    }

    // 找出間格值陣列 isovalues 在 [a,b] 之間的元素陣列
    vector<double>  find_subvalues( double a , double b , 
                                    const vector<double>& isovalues ) const {
        vector<double>  foo ;

        double  max0 = ( a > b ? a : b ) ;
        double  min0 = ( a > b ? b : a ) ;
        
        for ( int i = 0 ; i < isovalues.size() ; ++i ) {
            if ( min0 <= isovalues[i] && isovalues[i] <= max0 ) 
                foo.push_back(isovalues[i]) ;
        }

        return  foo ;

    }
    
    
  public :

    // 建構函式
    Triangular_Cell( const vector<Point>& p , const vector<double>& v ) 
        : Cell(p,v) {

        // 找出 deg
        deg = 1 ;
        while ( (deg+1)*(deg+2)/2 < no_pts ) ++deg ;

        // 根據 deg 設定使用的 shape function
        if ( shape_fn.size() == 0 ) setup_shape_function(deg) ;

    }

    // 回傳 no of edge
    int  no_edge() const { return  static_cast<int>(No_Edge) ; }

    // 回傳格子的總點數
    int  no_cell_pts() const { return (deg+1)*(deg+2)/2 ; }

    // 功能與上一樣，但設為 靜態函式
    static  int  no_cell_pts(int d) { return (d+1)*(d+2)/2 ; }
    

    // 回傳本格子的四個邊界點實際位置
    vector<Point>  cell_border_pts() const {
        return  vector<Point>(pts.begin(),pts.begin()+No_Edge) ;
    }
    
        
    // 將數值 [a,b] 區間設定 n 個等份， 輸出等數值限到檔案 fn 
    // Cell 的每邊切割成 m 個小等份，
    void  plot_isogram( ostream& outfile , double a , double b , int n , int m ,
                        bool  plot_cell_boundary ) const {

        int             i , j , k , s , t ;
        double          dx = 2./m ;
        vector<Point>   division_pts , border_pts ;
        vector<double>  iso_values(n+1) ;
        vector<double>  division_vals ;

        vector<Point>   isoline ;
        vector<double>  sub_isovals ;


        // 計算所要畫出的所有等值點
        double          dv = (b-a)/n ;
        for ( i = 0 ; i <= n ; ++i ) iso_values[i] = a + i*dv ;


        // 畫格子的邊界
        if ( plot_cell_boundary ) {
            border_pts = cell_border_pts() ;
            copy(border_pts.begin(),border_pts.end(),ostream_iterator<Point>(outfile,"\n")) ;
            outfile << border_pts[0] << "\n\n" ;
        }
        

        // 將格子每邊各切割成 m 等份，然後畫出每一個小區間範圍的等值線
        //
        for ( i = 0 ; i <= 2*(m-1) ; ++i ) {

            for ( j = (i%2?1:0) ; j <= 2*(m-1)-i ; j+=2 ) {

                // 找出圍住 (i,j) 小區間四個頂點的實際位置與其對應數值
                division_pts  = get_grid_points(i,j,m) ;
                division_vals = get_grid_val(i,j,m) ;

                // 找出每個小區間的 s 邊與剩餘的 t 邊的等值線的兩個端點
                for ( s = 0 ; s < No_Edge ; ++s ) {

                    // 找出在小區間的 s 邊上的實際可用的等值陣列
                    sub_isovals = find_subvalues(division_vals[s],
                                                 division_vals[(s+1)%4],
                                                 iso_values) ;

                    // 找出 s 到 t 邊的等值線
                    for ( t = (s+1)%4 ; t < No_Edge ; ++t ) {

                        for ( k = 0 ; k < sub_isovals.size() ; ++k ) {
                            isoline = division_isoline(division_pts,division_vals,
                                                       s,t,sub_isovals[k]) ;
                            if ( isoline.size() == 2 ) {
                                outfile << isoline[0] << "\n" << isoline[1] << "\n\n" ;
                            }
                        }

                        if ( s == No_Edge-1 ) break ;
                        
                    }
                }

            }
        }
        
    }


    ~Triangular_Cell(){}
    
};


vector<Triangular_Shape_Function*>  Triangular_Cell::shape_fn ;


class  Domain {
  private :
    
    string  ofile ;          // 輸出的資料檔
    double  vmin , vmax ;    // 設定最小與最大的數值
    int     val_division ;   // 數值切割的區間數
    int     cell_division ;  // cell 切割的區間數

    int     max_node ;       // 可能的最多點數

    mutable ofstream  outfile ;      // 設定輸出檔
    
    vector<Cell*>     cells ;

    // 讀取資料，井號開始的列為註解列
    void  read_data( char* filename ) {

        istringstream   istr ;
        string          line ;
        
        ifstream  infile(filename) ;   // 開啟資料檔

        // 讀取輸出檔名
        while ( getline(infile,line) ) {
            if ( line.size() > 0 ) {
                if ( line[0] == '#' ) continue ;
                istr.str(line) ;
                istr >> ofile ;
                outfile.open(ofile.c_str()) ;
                break ;
            }
        }
        istr.clear() ;

        // 讀取非格子的基本資料
        int     i = 0 ;
        double  tmp[5] ;
        while ( getline(infile,line) ) {
            if ( line.size() > 0 ) {
                if ( line[0] == '#' ) continue ;
                istr.str(line) ;
                while ( istr >> tmp[i] ) ++i ;
                istr.clear() ;
                if ( i == 5 ) break ;
            }
        }


        vmin = tmp[0] ;
        vmax = tmp[1] ;
        val_division  = static_cast<int>(tmp[2]) ;
        cell_division = static_cast<int>(tmp[3]) ;

        // 插分多項式的階數
        int  deg = static_cast<int>(tmp[4]) ;

        Point  p ;
        double v ;

        vector<Point>   pts ;
        vector<double>  vals ;

        int      cell_no ;             // 格子的數量
        int      no_pts ;              // 格子點的數量
        int      total_edge_no = 0 ;   // 總邊數
        
        cell_no = no_pts = 0 ;

        // 每次讀一列資料
        while ( getline(infile,line) ) {

            // 井號在第一個字元，則為空白行
            if ( line.size() != 0 && line[0] == '#' ) continue ;
            
            if ( line.find_first_not_of(" ") != string::npos ) {

                // 如果非空白字元，則讀入格子點資料
                istr.str(line) ;
                istr >> p >> v ;

                pts.push_back(p) ;
                vals.push_back(v) ;
                istr.clear() ;

                ++no_pts ;
                
            } else {

                if ( no_pts > 0 ) {

                    // 如果為空白字元且點數不為零，則產生格子物件

                    ++cell_no ;
                    
                    if ( no_pts == Rectangular_Cell::no_cell_pts(deg) ) {
                        cells.push_back( new Rectangular_Cell(pts,vals) ) ;
                    } else if ( no_pts == Triangular_Cell::no_cell_pts(deg) ) {
                        cells.push_back( new Triangular_Cell(pts,vals) ) ;
                    } else {
                        cout << "> 資料檔的第 " << cell_no
                             << " 筆格子點的點數資料錯誤" << endl ;
                        exit(1) ;
                    }

                    total_edge_no += cells.back()->no_edge() ;

                    no_pts = 0 ;
                    
                    // 如果為空白字元，則重設 pts 與 vals
                    pts.resize(0) ;
                    vals.resize(0) ;

                }

            }
            
        }

        // 最多可能的點數
        max_node = (cell_no+2) * total_edge_no ;
        
    }

    // 給某邊的兩端點的編號 i 與 j ，計算出邊的編號
    int  edge_no( int i , int j ) const {
        return i > j ? j*max_node+i : i*max_node+j ;
    }

    // 輸入邊的編號，回傳兩端點的編號
    pair<int,int>  edge_nodes( int i ) const {
        return make_pair<int,int>(i/max_node,i%max_node) ;
    }

    // 畫整個區域邊界線
    void  plot_domain_boundary() const {

        int             c , i , j , no ;

        map<Point,int>  pts ;      // 點 --> 編號
        map<int,int>    edges ;    // 邊的重複數

        vector<Point>  pts2(max_node) ;  // 編號 --> 點
        vector<Point>  bpts ;

        pair<int,int>   edge_node ;
        map<int,int>::iterator  iter ;
        
        int             a , b ;
        
        no = 1 ;
        for ( i = 0 ; i < cells.size() ; ++i ) {
            bpts = cells[i]->cell_border_pts() ;

            for ( j = 0 ; j < bpts.size() ; ++j ) {
                if ( pts[bpts[j]] == 0 ) {
                    pts[bpts[j]] = no ;
                    pts2[no] = bpts[j] ;
                    ++no ;
                }
            }

            for ( j = 0 ; j < cells[i]->no_edge() ; ++j ) {
                a = pts[bpts[j]] ;
                b = pts[bpts[(j+1)%cells[i]->no_edge()]] ;
                edges[edge_no(a,b)]++ ;
            }

        }

        outfile << "# 整個區域的邊界線\n" ;
        for ( iter = edges.begin() ; iter != edges.end() ; ++iter ) {
            if ( iter->second == 1 ) {
                edge_node = edge_nodes(iter->first) ;
                outfile << pts2[edge_node.first] << "\n" 
                        << pts2[edge_node.second] << "\n\n" ;
            }
        }

    }

    
  public :

    Domain( char* filename ) { read_data(filename) ; }
    
    // 以 gnuplot 格式輸出等值線
    void  plot_isogram() const {

        // 畫整個區域的邊界線
        plot_domain_boundary() ;
        
        // 是否要畫每個格子的邊界線
        bool  plot_cell_boundary = false ;

        // 畫每個格子的等值線
        for ( int i = 0 ; i < cells.size() ; ++i ) {
            outfile << "# 第 " << i << " 個格子等值線" << "\n" ;
            cells[i]->plot_isogram(outfile,vmin,vmax,val_division,cell_division,
                                   plot_cell_boundary) ;
        }
        
    }


    ~Domain() {
        for ( int i = 0 ; i < cells.size() ; ++i ) delete cells[i] ;
    }
    
            
};



int main() {

    Domain  foo("fem.dat") ;

    foo.plot_isogram() ;
    
    return 0 ;
}


/************************************************************************ 
   OUTPUT    OUTPUT    OUTPUT    OUTPUT    OUTPUT    OUTPUT 
 ************************************************************************ 

 <a href="./ss.jpeg">以下為使用 gnuplot 讀入 aa.dat 輸出檔所畫出來的等值線圖形</a>

 ************************************************************************/