|
說 明 |
應用問題區首頁 |
讀入三角或四邊形網格的資料與點的數值，程式可畫出等值線
由 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 ;
}
|
輸 出 |
|---|
以下為使用 gnuplot 讀入 aa.dat 輸出檔所畫出來的等值線圖形 
| |
|