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