www.gusucode.com > 岩体稳定性计算 单个裂隙生成 三维模型 应用于边坡 大坝等岩质结构源码程序 > main.m

    function FCG_3D 
%Open directory tosave combined data 
fprintf('Select Folder to save fracture data\n\n'); 
path_name = uigetdir(); 
 
while path_name == 0 
 fprintf('Folder not found \n\n Select Folder to save data\n\n'); 
 path_name = uigetdir(); 
end 
 
Nf = input('Number of fractures in fracture cluster: '); 
fprintf('\n'); 
while isempty(Nf) == 1 || rem(Nf,1) ~= 0 
 fprintf('need to specify a whole number\n\n'); 
 Nf = input('Number of fractures in fracture set: '); 
 fprintf('\n'); 
end 
 
Xmax = input('Horizontal length of domain (m): '); 
fprintf('\n'); 
while isempty(Xmax) == 1 
 fprintf('need to specify a number\n\n'); 
 Xmax = input('Horizontal length of domain (m): '); 
 fprintf('\n'); 
end 
 
Ymax = input('Width of domain into page (m): '); 
fprintf('\n'); 
while isempty(Ymax) == 1 
 fprintf('need to specify a number\n\n'); 
 Ymax = input('Width of domain into page (m): '); 
 fprintf('\n'); 
end 
 
Zmax = input('Vertical length of domain (m): '); 
fprintf('\n'); 
while isempty(Zmax) == 1 
 fprintf('need to specify a number\n\n'); 
 Zmax = input('Vertical length of domain (m): '); 
 fprintf('\n'); 
end 
 
xc = input('Horizontal location for each fracture (m): '); 
fprintf('\n'); 
while isempty(xc) == 1 || xc < 0 || xc > Xmax 
 fprintf('need to specify a number between 0 and Xmax \n\n'); 
 xc = input('Horizontal location for each fracture (m): '); 
 fprintf('\n'); 
end 
 
yc = input('Width location for each fracture (m): '); 
fprintf('\n'); 
while isempty(yc) == 1 || yc < 0 || yc > Ymax 
 fprintf('need to specify a number between 0 and Ymax \n\n'); 
 yc = input('Width location for each fracture (m): '); 
 fprintf('\n'); 
end 
 
 
zc = input('Vertical location for each fracture (m): '); 
fprintf('\n'); 
while isempty(zc) == 1 || zc < 0 || zc > Zmax 
 fprintf('need to specify a number between 0 and Zmax \n\n'); 
 zc = input('Vertical location for each fracture (m): '); 
 fprintf('\n'); 
end 
 
Npoints = input('Number of sides of fracture polygon: '); 
fprintf('\n'); 
while isempty(Npoints) == 1 || rem(Npoints,1) ~= 0 || Npoints < 3 
  fprintf('need to specify a whole number equal to or greater than 3 \n\n'); 
  Npoints = input('Number of sides of fracture polygon: '); 
  fprintf('\n'); 
end 
 
Mean_Size = input('Mean fracture size(m): '); 
fprintf('\n'); 
while isempty(Mean_Size) == 1 
 fprintf('need to specify a number\n\n'); 
 Mean_Size = input('Mean fracture size(m): '); 
 fprintf('\n'); 
end 
 
Dispersion_Size = input('Size standard deviation(m): '); 
fprintf('\n'); 
while isempty(Dispersion_Size) == 1 
 fprintf('need to specify a number\n\n'); 
 Dispersion_Size = input('Size standard deviation(m): '); 
 fprintf('\n'); 
end 
 
Shape_reg = input('fracture polygon shape to is  regular (type R) or general (type G): ','s'); 
fprintf('\n'); 
while size(Shape_reg,2)~= 1 || (Shape_reg ~= 'G' && Shape_reg ~= 'R') 
 fprintf('type R (have regular polygon shape) '); 
 fprintf('type G (have standard polygon shape)\n\n'); 
 Shape_reg = input('fracture polygon shape to is  regular (type R) or general (type G): ','s'); 
end 
 
Xpoints = nan(Nf,Npoints + 1); 
Ypoints = nan(Nf,Npoints + 1); 
Zpoints = nan(Nf,Npoints + 1); 
Dip = nan(Nf,1); 
Thetha = nan(Nf,1); 
x = nan(Nf,1); 
y = nan(Nf,1); 
z = nan(Nf,1); 
Size = nan(Nf,1); 
Size_h =  nan(Nf,1); 
frac_index_no = nan(Nf,1); 
normal_vectors = nan(Nf,3); 
plane_eqns = nan(Nf,4); 
pressure = nan(1,6); 
 
%Generation of fracture parameters 
  for i = 1:1:Nf 
  %Orientation parameters - generated via uniform distribution  
    Dip(i) = 2*pi*rand(1); 
    Thetha(i) = 2*pi*rand(1); 
    
  %Size - both sides of ellipse generated via lognormal distribution 
    mu = log((Mean_Size^2)/sqrt(Dispersion_Size^2+Mean_Size^2)); 
    sigma = sqrt(log((Dispersion_Size^2)/(Mean_Size^2)+1)); 
    Size(i) = 0.5*lognrnd(mu,sigma); 
     
    if Shape_reg == 'R' 
    Size_h(i) = Size(i); 
    else 
    Size_h(i) = 0.5*lognrnd(mu,sigma); 
    end 
     
  %Location parameters - generated via uniform distribution 
    x(i) = xc; 
    y(i) = yc; 
    z(i) = zc; 
     
  %fracture index no 
    frac_index_no(i) = i;  
  end 
   
%Save fracture data in text file 
fracture_data = [frac_index_no,Size,Size_h,Thetha,Dip,x,y,z]; 
Ind_Frac_Data_filename = fullfile(path_name,'Fracture Data.txt'); 
save(Ind_Frac_Data_filename,'fracture_data','-ascii');   
   
%Generate 2D fracture polygon and fracture plane based on fracture parameters 
  for i = 1:1:Nf 
    [normal_vectors(i,:),plane_eqns(i,:),Xpoints(i,:),Ypoints(i,:),Zpoints(i,:)] = Polygon3D(Thetha(i),Dip(i),Size(i),Size_h(i),x(i),y(i),z(i),Npoints);       
  end 
  
%Set colour for each fracture based on fracture_size 
[frac_colour] = frac_colour_size(Size,Size_h,Nf); 
 
%cut fractures to fit domain specified: 
 %find fractures that intersect with boundary 
  [~,frac_boundary_points] = int_frac_boundary(Xpoints,Ypoints,Zpoints,Thetha,Dip,x,y,z,Xmax,Ymax,Zmax,Nf,pressure,normal_vectors,plane_eqns); 
  
 %truncate fracture points that are not within the specified domain  
 [Xadj,Yadj,Zadj] = fracture_boundary_truncation(Xpoints,Ypoints,Zpoints,Thetha,Dip,x,y,z,Nf,Npoints,frac_boundary_points,Xmax,Ymax,Zmax); 
  
%Plot fractures after boundary truncation 
  for i = 1:1:Nf 
   X_plot = Xadj(i,:); Y_plot = Yadj(i,:); Z_plot = Zadj(i,:); 
   X_plot(isnan(X_plot) == 1) = []; 
   Y_plot(isnan(Y_plot) == 1) = []; 
   Z_plot(isnan(Z_plot) == 1) = []; 
    
   plot3(X_plot,Y_plot,Z_plot,'color',frac_colour(i,:)); 
   patch(X_plot,Y_plot,Z_plot,frac_colour(i,:)); 
   hold on; 
  end 
 
%plot boundary lines do distinctly show domain 
boundary_lines(Xmax,Ymax,Zmax); 
 
%plot fracture network within domain specified 
axis([0 Xmax 0 Ymax 0 Zmax]);  
 
end