GMS建模
===
建立XYZ檔
---
**1. 開QGIS將所選DEM圖層轉成.XYZ檔**

.xyz檔就是將DTM每個點位轉成X軸、Y軸、Z軸顯示


內政部20M DTM 使用的CRS(大地地理座標) 為 TWD97-121 (EPSG:3826)


**2. .XYZ檔打開後看到第一行(286XXX)為經度(X軸)
第二行(276XXX)為緯度(Y軸)
第三行為高程(Z軸)**

**3. 將.XYZ檔用Excel打開編輯，將X、Y軸從原點(0,0)開始重新編輯**

**4. 排序X軸順序
設置表格後 > 檢視 > 凍結表格 > 凍結首欄
常用 > 排序與篩選 > 由小至大**

**5. 記住有多少筆資料**
如果有一行是XYZ軸記得去掉

**6. 將排序完的資料另存新檔(.xyz)**

**7. 最後加上一些東西即可**

TIN
BEGT
ACTIVETIN
ID  1
TNAM "20mtest" (檔案名稱)
THORID "0"
MAT 1
GUID 978bf581-e7f2-4894-954e-eb7ba01da9b0 (識別碼)
VERT 18733(第五點中的資料筆數)
X Y Z數據...

TRI 0
ENDT



**P.S. 20m DEM XY以20為等距 6m DEM 以6為等距**

輸入檔案
---
**開啟GMS > 左上File > Open 打開.xyz檔 > 右鍵Tin Data裡的檔案 > Z高程建立資料庫**


**從上方工具列找Display Options > TIN Data > 勾選Vertices顯示DTM節點**



**TINs > Triangulate將各節點相連形成不規則三角形 (點→線)**



建立2D Mesh Data
---

**在上方 materals建立分層資料還有顏色(ex.崩積層,泥岩,砂岩...等)**


**TINs > Convert to > 2D Mesh 將不規則三角形節點轉成網格 (線→面)**

**最上層地層顏色為網格顏色**

設置鑽井資料
---

**左邊工作欄右鍵 New > Borehole 建立鑽井資料**


**要移至平面才能開始操作**

**建設所要得井位**

**編輯井位資料**
根據設置井位的XY軸，透過之前的XYZ檔的高程去假設出地層

**Display Option > Borehole Data 去調整井位顯示**



**全部井位確定建好 > 上方Borehole > Auto-Assign Horizons將各地層作編號分類**


建立3D Mesh Data
---
**至少要有==6個==地下水井且分散面的四周(不能太集中)才能建立3D Mesh**


**頂層使用Tin data，Constant自訂義模型底部高程多少**

**生成3維網格細節**
* Interpolation Method (使用的算法)
1. Inverse Distance Weighted (反距離權重法)
2. Natural Neighbor (自然對數法)
真實井位多用反距離，真實井位少用自然對數

* Nodel Function (節點內插法)
1. Constant (根據各點位置內插)
2. Gradient Plane (根據頂層坡度內插)
3. Quadratic (根據網格密度品質內插)
表面水平用 Constant ， 表面起伏用 Gradient

* Meshing Options (網格型式)
1. Prevent Pinchouts (只生成Pyramids & Hexahedra)
2. Minimum Element Thickness (網格最小高度)
3. Refine Elements (固定網格型態)
a. Refine All Elements (全部模型只使用 四面體網格Tetrahedral)
b. Refine Layers With Pyramid Elements (全部模型只使用 五面體網格Pyramid)

* Subdivide Material Layers (各地層分層)
1. Target Layer Thickness (每個網格固定高度)
2. Number of Layers (每個地層分幾層)


**跑出3D Mesh圖**


單純看建模情形可分各地層1層來看，要導出確定模型再細部分多層
不然一次分太多層對電腦記憶體造成大量負擔

模型處理&輸出
---
**右鍵3D Mesh裡檔案 > 輸出3D Mesh (mesh1.txt)**



**用Excel打開剛導出的檔案 (mesh1.txt)**
Excel有資料筆數限制，如果資料超過一百萬筆，以每一百萬為界線將資料分割後處理。




**1. 將E4T與ND分開個別處理**



**2. E4T刪除,存e4t (.csv)**


**3. ND刪除,存node(.csv)**


**4. 將<font color="#8F16C6">GMS2Flac3D</font> (Matlab Code)放進e4t(.csv)與node(.csv)同個資料夾並執行**

==**以下為全四面體轉換程式碼**== (有勾選Refine element)

function GMS2Flac3D
%從flac3d中取得資訊
tic
[xyz,elem]=read_GMS_files;
disp('GMS files read...')
toc
[k]=read(elem,xyz);
toc

save workspace
%diary('D:\FORMAT\lol.txt');diary on;

%% *******************************
%% *         FUNCTIONS           *
%% *******************************

function [xyz,elem]=read_GMS_files
    fprintf('node.csv：\n');
    type node.csv		
    xyz = csvread('node.csv');
    fprintf('e4t.csv：\n');
    type e4t.csv		
    elem = csvread('e4t.csv');
end

function [k]=read(elem,xyz)
    fid=fopen('GMS2FLAC.flac3d','a+');
    fprintf(fid,'* FLAC3D grid produced by FLAC3D\n');
    fprintf(fid,'* GRIDPOINTS\n');
    [m,~]=size(xyz);   
        for i=1:m
            nn=xyz(i,1);
		    xx=xyz(i,2);
		    yy=xyz(i,3);
		    zz=xyz(i,4);
            fprintf(fid,'G  %d   %e   %e   %e\n',nn,xx,yy,zz);
        end
    fprintf(fid,'* ZONES\n');
    [n,~]=size(elem);
        for i=1:n
            fN=elem(i,1);
            f0=elem(i,2);
            f1=elem(i,3);
            f2=elem(i,4);
            f3=elem(i,5);
            fprintf(fid,'Z T4   %d %d %d %d %d\n',fN,f0,f1,f2,f3);
            clear n
        end
    fprintf(fid,'* GROUPS\n');
    qq=6;  %NUMBER OF GROUPS
    [h,~]=size(elem);
        for i=1:qq
            G=i;
            fprintf(fid,'ZGROUP %d\n',G); 
            for j=1:h  
                k=[];
                GR=num2str(elem(j,6));
                bb=num2str(i);
                if GR==bb
                    w=(elem(j,1));  
                    k=(w);      
                    fprintf(fid,'%d\n',k);
                    j=j+1;
                else
                    disp('false')
                end
            end
        end
        clear h
    end
end

==**以下為四面、五面、六面體轉換程式**== (無勾選Refine element)
五面體與六面體網格處理方法同E4T

function GMS2Flac3D
%從flac3d中取得資訊
tic
[ND,e4t,e5p,e6w]=read_GMS_files;
disp('GMS files read...')
toc
[k]=read(e4t,e5p,e6w,ND);
toc
save workspace

%% *******************************
%% *         FUNCTIONS           *
%% *******************************

function [ND,e4t,e5p,e6w]=read_GMS_files
    fprintf('node.csv：\n');
    type node.csv
    ND = csvread('node.csv');
    fprintf('e4t.csv：\n');
    type e4t.csv
    e4t = csvread('e4t.csv');
    fprintf('e5p.csv：\n');
    type e5p.csv
    e5p = csvread('e5p.csv');
    fprintf('e6w.csv：\n');
    type e6w.csv
    e6w = csvread('e6w.csv');
end

function [k]=read(e4t,e5p,e6w,ND)
    fid=fopen('GMS2FLAC.flac3d','a+');
    fprintf(fid,'* FLAC3D grid produced by FLAC3D\n');
    fprintf(fid,'* GRIDPOINTS\n');

    [gridpoint,~]=size(ND);   
        for i=1:gridpoint
            nn=ND(i,1);
		    xx=ND(i,2);
		    yy=ND(i,3);
		    zz=ND(i,4);
            fprintf(fid,'G  %d   %e   %e   %e\n',nn,xx,yy,zz);
			clear gridpoint
        end

    fprintf(fid,'* ZONES\n');
    [tetrahedral,~]=size(e4t);
        for i=1:tetrahedral
            fN=e4t(i,1); %ID
            f0=e4t(i,2);
            f1=e4t(i,3);
            f2=e4t(i,4);
            f3=e4t(i,5);
            fprintf(fid,'Z T4   %d %d %d %d %d\n',fN,f0,f1,f2,f3);
            clear tetrahedral
        end

	[pyramid,~]=size(e5p);
        for i=1:pyramid
            fN=e5p(i,1); %ID
            f0=e5p(i,2);
            f1=e5p(i,3);
            f2=e5p(i,4);
            f3=e5p(i,5);
			f4=e5p(i,6);
            fprintf(fid,'Z P5   %d %d %d %d %d %d\n',fN,f0,f1,f3,f4,f2);
            clear pyramid
        end
		
	[wedge,~]=size(e6w);
        for i=1:wedge
            fN=e6w(i,1); %ID
            f0=e6w(i,2);
            f1=e6w(i,3);
            f2=e6w(i,4);
            f3=e6w(i,5);
			f4=e6w(i,6);
			f5=e6w(i,7);
            fprintf(fid,'Z W6   %d %d %d %d %d %d %d\n',fN,f0,f2,f3,f1,f5,f4);
            clear wedge
        end
		
    fprintf(fid,'* GROUPS\n');
    strata=4;  %地層個數

    [t4f,~]=size(e4t);
    [p5e,~]=size(e5p);
    [w6e,~]=size(e6w);
        for i=1:strata
            G=i;
            fprintf(fid,'ZGROUP %d\n',G); 
            for j=1:t4f  
                k=[];
                GR=num2str(e4t(j,6));
                bb=num2str(i);
                if GR==bb
                    w=(e4t(j,1));  
                    k=(w);      
                    fprintf(fid,'%d\n',k);
                    j=j+1;
                else
                    disp('false')
                end
            end
            for j=1:p5e  
                k=[];
                GR=num2str(e5p(j,7));
                bb=num2str(i);
                if GR==bb
                    w=(e5p(j,1));  
                    k=(w);      
                    fprintf(fid,'%d\n',k);
                    j=j+1;
                else
                    disp('false')
                end
            end
            for j=1:w6e  
                k=[];
                GR=num2str(e6w(j,8));
                bb=num2str(i);
                if GR==bb
                    w=(e6w(j,1));  
                    k=(w);      
                    fprintf(fid,'%d\n',k);
                    j=j+1;
                else
                    disp('false')
                end
            end                        
        end
    clear t4f
    clear p5e
    clear w6e
    end
end

**5. 跑完會生成<font color="#8F16C6">GMS2FLAC.flac3D</font>(如果要重跑要刪除後再跑)**

模型導入Flac3D
---

**將生成的GMS2FLAC.flac3d 用Flac3D 打開便能看見建立的模型**
左上File > Open > GMS2FLAC.flac3d



**打開Model視窗可以看到自己建立的模型**


**點選地層 > 工具列 Assign a Group... > 將所有地層做分類**



**工具列Assign Group names to face... > Ignore existing group names分東西南北面**

分面是要方便後面Flac3d做力學運算



**回到Plot視窗 > 右上+Zone可以看見分完地層的模型**




**==看到模型沒問題在下方輸入 model save 'ini_model' 將模型存檔==**

生成地下水位面
---

**1. 在borehole裡輸入地下水位高程**



**2. borehole(上面) > 選advanced > waterable... (幫他取名字water之類的)**



**3. 對Scatter右鍵 > interpolate to > active tin (轉成.tin檔,所以要在Tin Data裡面看)**


**4. 將地下水位.tin檔 > 生成2D Mesh**

**5. 2D Mesh > 生成 Cad圖檔(.dxf) > 將Cad圖檔存至運算的Flac3D資料夾內**

**6. 透過Flac3D程式碼便可在Water Table視窗看到地下水位面**

地下水位面(舊方法)
---
**4. 接續上面第3點到Tin Data裡將地層與地下水位資料輸出**


要輸出地層先點default再點surface輸出.tin檔


要輸出水位先點watertable再點surface輸出.tin檔


將輸出的地下水位面資料(surface.tin)用Matlab轉換成Flac3D可讀取的格式

clc
clear

X0_=0; %X軸原點
Y0_=0; %Y軸原點
BOT=50; %Z軸原點
xsize_=615; %X軸邊長
ysize_=355; %Y軸邊長
interval_=5; %節點間距(網格大小,DTM精度)

yline_=(ysize_/interval_)+1;
xline_=(xsize_/interval_)+1;
todata=yline_*xline_;
load('surface.tin');
XX=surface(:,1);
YY=surface(:,2);
TT=surface(:,3);
X1=X0_:interval_:xsize_;
Y1=Y0_:interval_:ysize_;
[X2,Y2]=meshgrid(X1,Y1);
Z1=griddata(XX,YY,TT,X2,Y2);
surf(X1,Y1,Z1)

%將文件轉換成Flac3D
for i=1:yline_
    no_(i,1:xline_)=i+100;
    BOTNO(i,1:xline_)=i+100;
    BOT_(i,1:xline_)=BOT;
end
DATA_(:,1)=flipud(rot90(reshape(rot90(no_), 1, todata)));
XD=flipud(rot90(X2));
ZD=flipud(rot90(Z1));
DATA_(:,2)=reshape(XD, 1, todata);
DATA_(:,3)=reshape(ZD, 1, todata);  
DATA_(isnan(DATA_(:,3)),3)=-0.5;
DATAB_(:,1)=flipud(rot90(reshape(rot90(BOTNO), 1, todata)));
DATAB_(:,2)=DATA_(:,2);
DATAB_(:,3)=reshape(BOT_, 1, todata);
ch(1:todata,1)=char('t');
ch(1:todata,2)=char('a');
ch(1:todata,3)=char('b');
ch(1:todata,4)=char('l');
ch(1:todata,5)=char('e');
ch(1:todata,6)=char(' ');
ftable1_=num2str(DATA_(:,1));
fx1_=num2str(DATA_(:,2));
fz1_=num2str(DATA_(:,3));
ftable2_=num2str(DATAB_(:,1));
fx2_=num2str(DATAB_(:,2));
fz2_=num2str(DATAB_(:,3));
temp={ch,ch(1:todata,6),ftable1_,ch(1:todata,6),fx1_,ch(1:todata,6),fz1_};
temp=cell2mat(temp);
dlmwrite('watertab.dat',';Tables of interpolated data from contour map','delimiter','');
dlmwrite('watertab.dat',temp,'-append', 'delimiter','');
dlmwrite('watertab.dat','     ','delimiter','','-append')