diff --git a/.gitignore b/.gitignore index b1d1224..4638bdc 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +model/ + # Created by https://www.toptal.com/developers/gitignore/api/matlab,microsoftoffice,latex # Edit at https://www.toptal.com/developers/gitignore?templates=matlab,microsoftoffice,latex diff --git a/.vscode/settings.json b/.vscode/settings.json index 1afc075..0f27d6a 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -6,9 +6,6 @@ "Krevor", "upscaling" ], - "cSpell.ignoreWords": [ - "mian" - ], "cSpell.enableFiletypes": [ "markdown" ], diff --git a/demo.m b/demo.m new file mode 100644 index 0000000..61ffd95 --- /dev/null +++ b/demo.m @@ -0,0 +1,162 @@ +%% StrataTrapper demonstration + +ppool = parpool; % start default parpool (optional) + +%% Input params + +params = gen_params(); +options = gen_options(); + +%% Subscaling demo + +subscale_demo(params, options); + +%% Grid geometry + +grid = grid_demo(); +is_permeable = rand(grid.cells.num,1) > 0.1; +mask = is_permeable(1:1000); % process only a subset of cells + +%% Run StrataTrapper + +enable_waitbar = true; + +strata_trapped = strata_trapper(grid, mask, params, options, enable_waitbar); + +%% Visualize saturation functions + +curves_plot(mask, strata_trapped, params); + +%% OGS inputs generation + +export_fut = parfeval(backgroundPool,@()ogs_export(G,mask,strata_trapped,'model/'), 0); % async call to function + +ogs_export_base(grid, strata_trapped, params); % base model parameters + +%% helpers + +function subscale_demo(params, options) + +options.subscale = 20 * meter(); + +[~, sub_porosity, sub_permeability, sub_entry_pressures, ~] = downscale(... + [400,400,0.1].*meter(),... + params, options ... + ); + +fig = figure('WindowStyle','docked','Name','Sub-scaling demonstration'); +tiles = tiledlayout(fig,'flow',TileSpacing='compact',Padding='tight'); + + function slice_plot(ax,data,title_str,units_str) + imagesc(ax,data',... + 'YData',linspace(0,400,size(data,2)), ... + 'XData',linspace(0,400,size(data,1))... + ); + + axis(ax,'equal'); + xlabel(ax,'x, m'); + ylabel(ax,'y, m'); + + ax.XLimitMethod="tight"; + ax.YLimitMethod="tight"; + + title(ax,title_str); + cb = colorbar(ax,"eastoutside"); + ylabel(cb,units_str); + end + +slice_plot(nexttile(tiles),squeeze(sub_porosity(2,:,:)),'Porosity',''); + +slice_plot(nexttile(tiles),squeeze(sub_entry_pressures(2,:,:)./kilo()),'Entry pressure','kPa'); + +slice_plot(nexttile(tiles),squeeze(sub_permeability(2,:,:)./(milli()*darcy())),'Permeability','mD'); + +end + +function [grid] = grid_demo() +grid.cartDims = [82, 123, 36]; +grid.cells.num = prod(grid.cartDims); +grid.cells.indexMap = 1:grid.cells.num; +grid.DX = 400 * meter() * ones(grid.cells.num,1); +grid.DY = grid.DX; +grid.DZ = 0.1 * meter() * ones(grid.cells.num,1); +end + +function curves_plot(mask, strata_trapped, params) +sub_data = @(data,mask,direction) squeeze(data(mask,direction,:)); + +fig_krw = figure('WindowStyle','docked','Name','krw'); +tiles_krw = tiledlayout(fig_krw,'flow',TileSpacing='compact',Padding='tight'); + +stat_plot(nexttile(tiles_krw),'Water relative permeability along X-axis','',strata_trapped.saturation,params.rel_perm.wat_func,sub_data(strata_trapped.rel_perm_wat,mask,1)); +stat_plot(nexttile(tiles_krw),'Water relative permeability along Y-axis','',strata_trapped.saturation,params.rel_perm.wat_func,sub_data(strata_trapped.rel_perm_wat,mask,2)); +stat_plot(nexttile(tiles_krw),'Water relative permeability along Z-axis','',strata_trapped.saturation,params.rel_perm.wat_func,sub_data(strata_trapped.rel_perm_wat,mask,3)); + +fig_krg = figure('WindowStyle','docked','Name','krg'); +tiles_krg = tiledlayout(fig_krg,'flow',TileSpacing='compact',Padding='tight'); + +stat_plot(nexttile(tiles_krg),'Gas relative permeability along X-axis','',strata_trapped.saturation,@(sw) params.rel_perm.gas_func(1-sw),sub_data(strata_trapped.rel_perm_gas,mask,1)); +stat_plot(nexttile(tiles_krg),'Gas relative permeability along Y-axis','',strata_trapped.saturation,@(sw) params.rel_perm.gas_func(1-sw),sub_data(strata_trapped.rel_perm_gas,mask,2)); +stat_plot(nexttile(tiles_krg),'Gas relative permeability along Z-axis','',strata_trapped.saturation,@(sw) params.rel_perm.gas_func(1-sw),sub_data(strata_trapped.rel_perm_gas,mask,3)); + +fig_pc = figure('WindowStyle','docked','Name','pc'); +stat_plot(axes(fig_pc),'Capillary pressure','bar',strata_trapped.saturation,@(sw) params.capil.pres_func(sw,params.capil.pres_entry)./barsa(),strata_trapped.capillary_pressure(mask,:)./barsa()); + +end + +function ax = stat_plot(ax, name, y_label, x_data, base_func, data) +arguments + ax + name char + y_label char + x_data (1,:) double + base_func + data (:,:) double +end + +parallelcoords(ax,data,'Quantile',0.01,'XData',x_data); + +if ~isempty(base_func) + hold(ax,'on'); + plot(x_data,base_func(x_data),'-r'); + hold(ax,'off'); +end + +title(ax,name); + +xlabel(ax,'Wetting phase saturation'); +xlim(ax,[0,1]); +ax.XTickMode='auto'; +ax.XTickLabelMode='auto'; +ax.XLimitMethod="tickaligned"; + +% ylim(ax,[0,max(data,[],'all')]); +ylabel(ax,y_label); +ax.YLimitMethod="tight"; + +legend(ax,{'Median','Quantiles 0.01 and 0.99','','Intrinsic curve'},'Location','best'); +end + +function ogs_export_base(G, strata_trapped, params) +perm_base = params.perm_corr(strata_trapped.porosity)./milli()./darcy(); +perm_base(strata_trapped.porosity == 0) = 0; + +perm_base_mapping = zeros(G.cells.num,1); +perm_base_mapping(G.cells.indexMap) = perm_base; + +write_curve_nums("model/base-PERMX.grdecl","PERMX",perm_base_mapping,0,0); +write_curve_nums("model/base-PERMY.grdecl","PERMY",perm_base_mapping,0,0); +write_curve_nums("model/base-PERMZ.grdecl","PERMZ",perm_base_mapping,0,0); + +swfn_base = [strata_trapped.saturation; ... + params.rel_perm.wat_func(strata_trapped.saturation); ... + params.capil.pres_func(strata_trapped.saturation,params.capil.pres_entry)./barsa()]; + +sg = unique([1-strata_trapped.saturation,1]); + +sgfn_base = [sg; params.rel_perm.gas_func(sg); zeros(size(sg))]; + +format_spec ='%f\t%f\t%f\n'; +swfn_base_str = sprintf(format_spec,swfn_base) +sgfn_base_str = sprintf(format_spec,sgfn_base) +end diff --git a/ogs_export.m b/ogs_export.m new file mode 100644 index 0000000..e235150 --- /dev/null +++ b/ogs_export.m @@ -0,0 +1,19 @@ +function [] = ogs_export(G, mask, strata_trapped, output_prefix) +arguments + G (1,1) struct + mask (:,1) logical + strata_trapped (1,1) struct + output_prefix char = 'model/' +end + +write_mappings(output_prefix,G,strata_trapped.permeability ./ milli ./ darcy(), strata_trapped.porosity); + +generate_sfn(... + mask, ... + strata_trapped.saturation, ... + strata_trapped.capillary_pressure ./ barsa(), ... + strata_trapped.rel_perm_wat, ... + strata_trapped.rel_perm_gas, ... + output_prefix, '.dat'); +end + diff --git a/src/brooks_corey_inv.m b/src/brooks_corey_inv.m new file mode 100644 index 0000000..5f9b211 --- /dev/null +++ b/src/brooks_corey_inv.m @@ -0,0 +1,3 @@ +function sw = brooks_corey_inv(pc,sw_resid,p_entry,lambda) +sw = (p_entry./pc).^lambda .* (1-sw_resid) + sw_resid; +end \ No newline at end of file diff --git a/src/calc_entry_pressure.m b/src/calc_entry_pressure.m new file mode 100644 index 0000000..ad44710 --- /dev/null +++ b/src/calc_entry_pressure.m @@ -0,0 +1,11 @@ +function entry_pressure = calc_entry_pressure(poro, perm, j_entry, contact_angle, surface_tension) +arguments + poro double + perm double + j_entry (1,1) double + contact_angle (1,1) double + surface_tension (1,1) double +end + +entry_pressure = j_entry * cos(contact_angle) * surface_tension * sqrt(poro./perm); +end \ No newline at end of file diff --git a/Functions/check_axis_connectivity.m b/src/check_axis_connectivity.m similarity index 97% rename from Functions/check_axis_connectivity.m rename to src/check_axis_connectivity.m index e9165d2..1269b9d 100644 --- a/Functions/check_axis_connectivity.m +++ b/src/check_axis_connectivity.m @@ -1,2101 +1,2101 @@ -function [connected, connected_mat] = check_axis_connectivity(kr_phase_mat, Nx, Nz,Ny,tol, direction, end_val) - - -connected_mat = zeros(Nz, Nx, Ny); - -if (direction == 1) - - % disp("Direction is 1"); - - %Direction 1 = 2D surface in x-y, moving down through z - - changing = 1; - connected = 0; - connected_mat(:,:,1) = 1; - - - while (changing == 1) - n_changed = 0; - - for ii = 1:Nz - - if (connected == 0) - for jj = 1:Nx - for kk = 1:Ny - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %Scenarios where 3 faces are blocked (8 cases when this can - %happen) - - if (ii == 1) && (jj == 1) && (kk == 1) %starting cube, 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - - elseif (ii == Nz) && (jj == Nx) && (kk == Ny) %end cube, 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - - elseif (ii == Nz) && (jj == 1) && (kk == 1) % (N, 1, 1), 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - elseif (ii == 1) && (jj == 1) && (kk == Ny) %(1,1,N) , 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - elseif (ii == 1) && (jj == Nx) && (kk == 1) %(1,N,1) , 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - - elseif (ii == Nz) && (jj == Nx) && (kk == 1) %(N,N,1) , 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - elseif (ii == 1) && (jj == Nx) && (kk == Ny) %(1,N,N), 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - elseif (ii == Nz) && (jj == 1) && (kk == Ny) %(N,1,N), 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %Scenario when 1 side is blocked (moving along 1 face, in the middle/centre box, no edges), 6 cases - - - - elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == 1) %if on jj boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %jj boundary - - check_pos(1,1) = ii + 1; - check_pos(2,1) = ii - 1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj+1; %only option for jj to increase is to go up - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - - elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == 1) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %ii boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii+1; %only option for ii to increase is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == 1) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %kk boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - check_pos(5,1) = ii; %only option for ii to increase is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - check_pos(5,3) = kk+1; %only way is for kk to increase - - elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == Nx) %if on jj boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %jj boundary - - check_pos(1,1) = ii + 1; - check_pos(2,1) = ii - 1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj-1; %only option for jj to decrease is to go up - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - - elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == Nz) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %ii boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii-1; %only option for ii to decrease is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == Ny) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %kk boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - check_pos(5,1) = ii; %only option for ii to increase is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - check_pos(5,3) = kk-1; %only way is for kk to increase - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %Now: 2 sides blocked (edges of the faces), 12 cases - - %first 4: k is in middle - elseif (ii == 1) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - elseif (ii == Nz) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - elseif (ii == Nz) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - - elseif (ii == 1) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - %next 4: i is in middle - - elseif (kk == 1) && (jj == 1) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (jj == Nx) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (jj == 1) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == 1) && (jj == Nx) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - %next 4 cases: j is in middle - - elseif (kk == 1) && (ii == 1) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii+1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (ii == Nz) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii-1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (ii == 1) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii+1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == 1) && (ii == Nz) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii-1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %centre cubes - no sides are blocked, only need 1 case - - - elseif ((ii > 1) && (ii < Nz) && (jj >1) && (jj < Nx) && (kk > 1) && (kk < Ny)) - - - n_check = 6; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - check_pos(5,1) = ii; - check_pos(6,1) = ii; - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - check_pos(6,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - check_pos(5,3) = kk+1; - check_pos(6,3) = kk-1; - - else - - - - disp("Error- cant find a situation") - - disp(ii) - disp(jj) - disp(kk) - - disp(Nz) - disp(Nx) - disp(Ny) - - - end - - %check if values are outside of our domain - % disp("Checking for errors") - - for i = 1:n_check - if (check_pos(i,2) > Nx) - - disp("i index exceeds!") - %disp(check_pos(:,:)) - % disp(ii) - % disp(jj) - % disp(kk) - % disp(Nx) - pause - end - - if (check_pos(i,1) > Nz) - disp("j index exceeds!") - % disp(ii) - % disp(jj) - % disp(kk) - - pause - end - - - if (check_pos(i,3) > Ny) - disp("k index exceeds!") - %disp(ii) - %disp(jj) - %disp(kk) - - pause - end - end - - %disp("movingon") - for i = 1:n_check - - %|| this means 'or' - % we are checking if we are at a boundary - - if (connected_mat(check_pos(i,1), check_pos(i,2), check_pos(i,3)) > 0) && (kr_phase_mat(ii,jj,kk) > tol) && (connected_mat(ii,jj,kk) == 0) - - connected_mat(ii,jj,kk) = 1; - n_changed = n_changed + 1; - end - - end - - - if (connected_mat(ii,jj,kk) == 1) && (kk == end_val) - connected = 1; - changing = 0; - break - else - connected = 0; - end - - - - - - end - - end - end - end - - if (n_changed == 0) - changing = 0; - end - - - end - - %%%%%%%%%%%%%%%%%""""""""""""""""""" DIRECTION 2 %%%%%%%%%%%%%%%%%"""""""""""""""""""%%%%%%%%%%%%%%%%%"""""""""""""""""""%%%%%%%%%%%%%%%%%""""""""""""""""""" - -elseif (direction == 2) - - % disp("Direction is 2"); - - %Direction 1 = 2D surface in x-z, moving down through y - - changing = 1; - connected = 0; - connected_mat(:,1,:) = 1; - - - while (changing == 1) - n_changed = 0; - - for kk = 1:Ny - - if (connected == 0) - for ii = 1:Nz - for jj = 1:Nx - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %Scenarios where 3 faces are blocked (8 cases when this can - %happen) - - if (ii == 1) && (jj == 1) && (kk == 1) %starting cube, 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - - elseif (ii == Nz) && (jj == Nx) && (kk == Ny) %end cube, 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - - elseif (ii == Nz) && (jj == 1) && (kk == 1) % (N, 1, 1), 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - elseif (ii == 1) && (jj == 1) && (kk == Ny) %(1,1,N) , 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - elseif (ii == 1) && (jj == Nx) && (kk == 1) %(1,N,1) , 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - - elseif (ii == Nz) && (jj == Nx) && (kk == 1) %(N,N,1) , 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - elseif (ii == 1) && (jj == Nx) && (kk == Ny) %(1,N,N), 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - elseif (ii == Nz) && (jj == 1) && (kk == Ny) %(N,1,N), 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %Scenario when 1 side is blocked (moving along 1 face, in the middle/centre box, no edges), 6 cases - - - - elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == 1) %if on jj boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %jj boundary - - check_pos(1,1) = ii + 1; - check_pos(2,1) = ii - 1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj+1; %only option for jj to increase is to go up - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - - elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == 1) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %ii boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii+1; %only option for ii to increase is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == 1) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %kk boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - check_pos(5,1) = ii; %only option for ii to increase is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - check_pos(5,3) = kk+1; %only way is for kk to increase - - elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == Nx) %if on jj boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %jj boundary - - check_pos(1,1) = ii + 1; - check_pos(2,1) = ii - 1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj-1; %only option for jj to decrease is to go up - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - - elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == Nz) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %ii boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii-1; %only option for ii to decrease is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == Ny) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %kk boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - check_pos(5,1) = ii; %only option for ii to increase is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - check_pos(5,3) = kk-1; %only way is for kk to increase - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %Now: 2 sides blocked (edges of the faces), 12 cases - - %first 4: k is in middle - elseif (ii == 1) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - elseif (ii == Nz) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - elseif (ii == Nz) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - - elseif (ii == 1) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - %next 4: i is in middle - - elseif (kk == 1) && (jj == 1) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (jj == Nx) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (jj == 1) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == 1) && (jj == Nx) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - %next 4 cases: j is in middle - - elseif (kk == 1) && (ii == 1) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii+1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (ii == Nz) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii-1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (ii == 1) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii+1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == 1) && (ii == Nz) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii-1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %centre cubes - no sides are blocked, only need 1 case - - - elseif ((ii > 1) && (ii < Nz) && (jj >1) && (jj < Nx) && (kk > 1) && (kk < Ny)) - - - n_check = 6; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - check_pos(5,1) = ii; - check_pos(6,1) = ii; - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - check_pos(6,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - check_pos(5,3) = kk+1; - check_pos(6,3) = kk-1; - - - - else - - - - disp("Error- cant find a situation") - - disp(ii) - disp(jj) - disp(kk) - - disp(Nz) - disp(Nx) - disp(Ny) - - - end - - %check if values are outside of our domain - - % disp("checking for errors - dir 2"); - - for i = 1:n_check - if (check_pos(i,2) > Nx) - disp("i out of range - dir 2") - pause - end - - if (check_pos(i,1) > Nz) - disp("j out of range - dir 2") - pause - end - - - if (check_pos(i,3) > Ny) - disp("k out of range - dir 3") - pause - end - end - - - for i = 1:n_check - - %|| this means 'or' - % we are checking if we are at a boundary - - if (connected_mat(check_pos(i,1), check_pos(i,2), check_pos(i,3)) > 0) && (kr_phase_mat(ii,jj,kk) > tol) && (connected_mat(ii,jj,kk) == 0) - - connected_mat(ii,jj,kk) = 1; - n_changed = n_changed + 1; - end - - end - - - if (connected_mat(ii,jj,kk) == 1) && (jj == end_val) - connected = 1; - changing = 0; - break - else - connected = 0; - end - - - - - - end - - end - end - end - - if (n_changed == 0) - changing = 0; - end - - - end - %%%%%%%%%%%%%%%%%""""""""""""""""""" DIRECTION 2 %%%%%%%%%%%%%%%%%"""""""""""""""""""%%%%%%%%%%%%%%%%%"""""""""""""""""""%%%%%%%%%%%%%%%%%""""""""""""""""""" - -elseif (direction == 3) - - %disp("Direction is 3"); - - %Direction 1 = 2D surface in y-z, moving down through x - - changing = 1; - connected = 0; - connected_mat(1,:,:) = 1; - - - while (changing == 1) - n_changed = 0; - - for jj = 1:Nx - - if (connected == 0) - for kk = 1:Ny - for ii = 1:Nz - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %Scenarios where 3 faces are blocked (8 cases when this can - %happen) - - if (ii == 1) && (jj == 1) && (kk == 1) %starting cube, 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - - elseif (ii == Nz) && (jj == Nx) && (kk == Ny) %end cube, 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - - elseif (ii == Nz) && (jj == 1) && (kk == 1) % (N, 1, 1), 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - elseif (ii == 1) && (jj == 1) && (kk == Ny) %(1,1,N) , 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - elseif (ii == 1) && (jj == Nx) && (kk == 1) %(1,N,1) , 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - - elseif (ii == Nz) && (jj == Nx) && (kk == 1) %(N,N,1) , 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - - elseif (ii == 1) && (jj == Nx) && (kk == Ny) %(1,N,N), 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - elseif (ii == Nz) && (jj == 1) && (kk == Ny) %(N,1,N), 3 faces blocked - - n_check = 3; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk-1; - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %Scenario when 1 side is blocked (moving along 1 face, in the middle/centre box, no edges), 6 cases - - - - elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == 1) %if on jj boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %jj boundary - - check_pos(1,1) = ii + 1; - check_pos(2,1) = ii - 1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj+1; %only option for jj to increase is to go up - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - - elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == 1) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %ii boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii+1; %only option for ii to increase is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == 1) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %kk boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - check_pos(5,1) = ii; %only option for ii to increase is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - check_pos(5,3) = kk+1; %only way is for kk to increase - - elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == Nx) %if on jj boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %jj boundary - - check_pos(1,1) = ii + 1; - check_pos(2,1) = ii - 1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj-1; %only option for jj to decrease is to go up - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - - elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == Nz) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %ii boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - check_pos(5,1) = ii-1; %only option for ii to decrease is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - check_pos(5,3) = kk; - - elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == Ny) %if on ii boundary, one face blocked off - - n_check = 5; - check_pos = zeros(n_check, 5); - - %kk boundary - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - check_pos(5,1) = ii; %only option for ii to increase is to go up - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - check_pos(5,3) = kk-1; %only way is for kk to increase - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %Now: 2 sides blocked (edges of the faces), 12 cases - - %first 4: k is in middle - elseif (ii == 1) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - elseif (ii == Nz) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - elseif (ii == Nz) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii-1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - - elseif (ii == 1) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off - - n_check = 4; - check_pos = zeros(n_check, 4); - - %kk boundary - - check_pos(1,1) = ii+1; - check_pos(2,1) = ii; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk+1; - check_pos(4,3) = kk-1; - - %next 4: i is in middle - - elseif (kk == 1) && (jj == 1) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (jj == Nx) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (jj == 1) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj+1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == 1) && (jj == Nx) && (ii > 1) && (ii < Nz) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - - check_pos(1,2) = jj; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - %next 4 cases: j is in middle - - elseif (kk == 1) && (ii == 1) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii+1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (ii == Nz) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii-1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == Ny) && (ii == 1) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii+1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk-1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - elseif (kk == 1) && (ii == Nz) && (jj > 1) && (jj < Nx) - - n_check = 4; - check_pos = zeros(n_check, 4); - - - check_pos(1,1) = ii; - check_pos(2,1) = ii-1; - check_pos(3,1) = ii; - check_pos(4,1) = ii; - - check_pos(1,2) = jj; - check_pos(2,2) = jj; - check_pos(3,2) = jj+1; - check_pos(4,2) = jj-1; - - check_pos(1,3) = kk+1; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %centre cubes - no sides are blocked, only need 1 case - - - elseif ((ii > 1) && (ii < Nz) && (jj >1) && (jj < Nx) && (kk > 1) && (kk < Ny)) - - - n_check = 6; - check_pos = zeros(n_check, 3); - - check_pos(1,1) = ii; - check_pos(2,1) = ii; - check_pos(3,1) = ii+1; - check_pos(4,1) = ii-1; - check_pos(5,1) = ii; - check_pos(6,1) = ii; - - check_pos(1,2) = jj+1; - check_pos(2,2) = jj-1; - check_pos(3,2) = jj; - check_pos(4,2) = jj; - check_pos(5,2) = jj; - check_pos(6,2) = jj; - - check_pos(1,3) = kk; - check_pos(2,3) = kk; - check_pos(3,3) = kk; - check_pos(4,3) = kk; - check_pos(5,3) = kk+1; - check_pos(6,3) = kk-1; - else - - disp("Error- cant find a situation") - - disp(ii) - disp(jj) - disp(kk) - - disp(Nz) - disp(Nx) - disp(Ny) - - end - - %check if values are outside of our domain - - %disp("checking for errors") - - for i = 1:n_check - if (check_pos(i,2) > Nx) - disp("i out of range") - - pause - end - - if (check_pos(i,1) > Nz) - disp("j out of range") - - pause - end - - - if (check_pos(i,3) > Ny) - disp("k out of range") - - pause - end - end - - - for i = 1:n_check - - %|| this means 'or' - % we are checking if we are at a boundary - - if (connected_mat(check_pos(i,1), check_pos(i,2), check_pos(i,3)) > 0) && (kr_phase_mat(ii,jj,kk) > tol) && (connected_mat(ii,jj,kk) == 0) - - connected_mat(ii,jj,kk) = 1; - n_changed = n_changed + 1; - end - - end - - - if (connected_mat(ii,jj,kk) == 1) && (ii == end_val) - connected = 1; - changing = 0; - break - else - connected = 0; - end - end - end - end - end - - if (n_changed == 0) - changing = 0; - end - - end -end -end +function [connected, connected_mat] = check_axis_connectivity(kr_phase_mat, Nx, Nz,Ny,tol, direction, end_val) + + +connected_mat = zeros(Nz, Nx, Ny); + +if (direction == 1) + + % disp("Direction is 1"); + + %Direction 1 = 2D surface in x-y, moving down through z + + changing = 1; + connected = 0; + connected_mat(:,:,1) = 1; + + + while (changing == 1) + n_changed = 0; + + for ii = 1:Nz + + if (connected == 0) + for jj = 1:Nx + for kk = 1:Ny + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %Scenarios where 3 faces are blocked (8 cases when this can + %happen) + + if (ii == 1) && (jj == 1) && (kk == 1) %starting cube, 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + + elseif (ii == Nz) && (jj == Nx) && (kk == Ny) %end cube, 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + + elseif (ii == Nz) && (jj == 1) && (kk == 1) % (N, 1, 1), 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + elseif (ii == 1) && (jj == 1) && (kk == Ny) %(1,1,N) , 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + elseif (ii == 1) && (jj == Nx) && (kk == 1) %(1,N,1) , 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + + elseif (ii == Nz) && (jj == Nx) && (kk == 1) %(N,N,1) , 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + elseif (ii == 1) && (jj == Nx) && (kk == Ny) %(1,N,N), 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + elseif (ii == Nz) && (jj == 1) && (kk == Ny) %(N,1,N), 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %Scenario when 1 side is blocked (moving along 1 face, in the middle/centre box, no edges), 6 cases + + + + elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == 1) %if on jj boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %jj boundary + + check_pos(1,1) = ii + 1; + check_pos(2,1) = ii - 1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj+1; %only option for jj to increase is to go up + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + + elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == 1) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %ii boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii+1; %only option for ii to increase is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == 1) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %kk boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + check_pos(5,1) = ii; %only option for ii to increase is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + check_pos(5,3) = kk+1; %only way is for kk to increase + + elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == Nx) %if on jj boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %jj boundary + + check_pos(1,1) = ii + 1; + check_pos(2,1) = ii - 1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj-1; %only option for jj to decrease is to go up + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + + elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == Nz) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %ii boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii-1; %only option for ii to decrease is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == Ny) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %kk boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + check_pos(5,1) = ii; %only option for ii to increase is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + check_pos(5,3) = kk-1; %only way is for kk to increase + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %Now: 2 sides blocked (edges of the faces), 12 cases + + %first 4: k is in middle + elseif (ii == 1) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + elseif (ii == Nz) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + elseif (ii == Nz) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + + elseif (ii == 1) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + %next 4: i is in middle + + elseif (kk == 1) && (jj == 1) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (jj == Nx) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (jj == 1) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == 1) && (jj == Nx) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + %next 4 cases: j is in middle + + elseif (kk == 1) && (ii == 1) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii+1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (ii == Nz) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii-1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (ii == 1) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii+1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == 1) && (ii == Nz) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii-1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %centre cubes - no sides are blocked, only need 1 case + + + elseif ((ii > 1) && (ii < Nz) && (jj >1) && (jj < Nx) && (kk > 1) && (kk < Ny)) + + + n_check = 6; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + check_pos(5,1) = ii; + check_pos(6,1) = ii; + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + check_pos(6,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + check_pos(5,3) = kk+1; + check_pos(6,3) = kk-1; + + else + + + + disp("Error- cant find a situation") + + disp(ii) + disp(jj) + disp(kk) + + disp(Nz) + disp(Nx) + disp(Ny) + + + end + + %check if values are outside of our domain + % disp("Checking for errors") + + for i = 1:n_check + if (check_pos(i,2) > Nx) + + disp("i index exceeds!") + %disp(check_pos(:,:)) + % disp(ii) + % disp(jj) + % disp(kk) + % disp(Nx) + pause + end + + if (check_pos(i,1) > Nz) + disp("j index exceeds!") + % disp(ii) + % disp(jj) + % disp(kk) + + pause + end + + + if (check_pos(i,3) > Ny) + disp("k index exceeds!") + %disp(ii) + %disp(jj) + %disp(kk) + + pause + end + end + + %disp("movingon") + for i = 1:n_check + + %|| this means 'or' + % we are checking if we are at a boundary + + if (connected_mat(check_pos(i,1), check_pos(i,2), check_pos(i,3)) > 0) && (kr_phase_mat(ii,jj,kk) > tol) && (connected_mat(ii,jj,kk) == 0) + + connected_mat(ii,jj,kk) = 1; + n_changed = n_changed + 1; + end + + end + + + if (connected_mat(ii,jj,kk) == 1) && (kk == end_val) + connected = 1; + changing = 0; + break + else + connected = 0; + end + + + + + + end + + end + end + end + + if (n_changed == 0) + changing = 0; + end + + + end + + %%%%%%%%%%%%%%%%%""""""""""""""""""" DIRECTION 2 %%%%%%%%%%%%%%%%%"""""""""""""""""""%%%%%%%%%%%%%%%%%"""""""""""""""""""%%%%%%%%%%%%%%%%%""""""""""""""""""" + +elseif (direction == 2) + + % disp("Direction is 2"); + + %Direction 1 = 2D surface in x-z, moving down through y + + changing = 1; + connected = 0; + connected_mat(:,1,:) = 1; + + + while (changing == 1) + n_changed = 0; + + for kk = 1:Ny + + if (connected == 0) + for ii = 1:Nz + for jj = 1:Nx + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %Scenarios where 3 faces are blocked (8 cases when this can + %happen) + + if (ii == 1) && (jj == 1) && (kk == 1) %starting cube, 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + + elseif (ii == Nz) && (jj == Nx) && (kk == Ny) %end cube, 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + + elseif (ii == Nz) && (jj == 1) && (kk == 1) % (N, 1, 1), 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + elseif (ii == 1) && (jj == 1) && (kk == Ny) %(1,1,N) , 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + elseif (ii == 1) && (jj == Nx) && (kk == 1) %(1,N,1) , 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + + elseif (ii == Nz) && (jj == Nx) && (kk == 1) %(N,N,1) , 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + elseif (ii == 1) && (jj == Nx) && (kk == Ny) %(1,N,N), 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + elseif (ii == Nz) && (jj == 1) && (kk == Ny) %(N,1,N), 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %Scenario when 1 side is blocked (moving along 1 face, in the middle/centre box, no edges), 6 cases + + + + elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == 1) %if on jj boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %jj boundary + + check_pos(1,1) = ii + 1; + check_pos(2,1) = ii - 1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj+1; %only option for jj to increase is to go up + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + + elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == 1) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %ii boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii+1; %only option for ii to increase is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == 1) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %kk boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + check_pos(5,1) = ii; %only option for ii to increase is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + check_pos(5,3) = kk+1; %only way is for kk to increase + + elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == Nx) %if on jj boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %jj boundary + + check_pos(1,1) = ii + 1; + check_pos(2,1) = ii - 1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj-1; %only option for jj to decrease is to go up + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + + elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == Nz) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %ii boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii-1; %only option for ii to decrease is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == Ny) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %kk boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + check_pos(5,1) = ii; %only option for ii to increase is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + check_pos(5,3) = kk-1; %only way is for kk to increase + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %Now: 2 sides blocked (edges of the faces), 12 cases + + %first 4: k is in middle + elseif (ii == 1) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + elseif (ii == Nz) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + elseif (ii == Nz) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + + elseif (ii == 1) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + %next 4: i is in middle + + elseif (kk == 1) && (jj == 1) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (jj == Nx) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (jj == 1) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == 1) && (jj == Nx) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + %next 4 cases: j is in middle + + elseif (kk == 1) && (ii == 1) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii+1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (ii == Nz) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii-1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (ii == 1) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii+1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == 1) && (ii == Nz) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii-1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %centre cubes - no sides are blocked, only need 1 case + + + elseif ((ii > 1) && (ii < Nz) && (jj >1) && (jj < Nx) && (kk > 1) && (kk < Ny)) + + + n_check = 6; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + check_pos(5,1) = ii; + check_pos(6,1) = ii; + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + check_pos(6,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + check_pos(5,3) = kk+1; + check_pos(6,3) = kk-1; + + + + else + + + + disp("Error- cant find a situation") + + disp(ii) + disp(jj) + disp(kk) + + disp(Nz) + disp(Nx) + disp(Ny) + + + end + + %check if values are outside of our domain + + % disp("checking for errors - dir 2"); + + for i = 1:n_check + if (check_pos(i,2) > Nx) + disp("i out of range - dir 2") + pause + end + + if (check_pos(i,1) > Nz) + disp("j out of range - dir 2") + pause + end + + + if (check_pos(i,3) > Ny) + disp("k out of range - dir 3") + pause + end + end + + + for i = 1:n_check + + %|| this means 'or' + % we are checking if we are at a boundary + + if (connected_mat(check_pos(i,1), check_pos(i,2), check_pos(i,3)) > 0) && (kr_phase_mat(ii,jj,kk) > tol) && (connected_mat(ii,jj,kk) == 0) + + connected_mat(ii,jj,kk) = 1; + n_changed = n_changed + 1; + end + + end + + + if (connected_mat(ii,jj,kk) == 1) && (jj == end_val) + connected = 1; + changing = 0; + break + else + connected = 0; + end + + + + + + end + + end + end + end + + if (n_changed == 0) + changing = 0; + end + + + end + %%%%%%%%%%%%%%%%%""""""""""""""""""" DIRECTION 2 %%%%%%%%%%%%%%%%%"""""""""""""""""""%%%%%%%%%%%%%%%%%"""""""""""""""""""%%%%%%%%%%%%%%%%%""""""""""""""""""" + +elseif (direction == 3) + + %disp("Direction is 3"); + + %Direction 1 = 2D surface in y-z, moving down through x + + changing = 1; + connected = 0; + connected_mat(1,:,:) = 1; + + + while (changing == 1) + n_changed = 0; + + for jj = 1:Nx + + if (connected == 0) + for kk = 1:Ny + for ii = 1:Nz + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %Scenarios where 3 faces are blocked (8 cases when this can + %happen) + + if (ii == 1) && (jj == 1) && (kk == 1) %starting cube, 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + + elseif (ii == Nz) && (jj == Nx) && (kk == Ny) %end cube, 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + + elseif (ii == Nz) && (jj == 1) && (kk == 1) % (N, 1, 1), 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + elseif (ii == 1) && (jj == 1) && (kk == Ny) %(1,1,N) , 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + elseif (ii == 1) && (jj == Nx) && (kk == 1) %(1,N,1) , 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + + elseif (ii == Nz) && (jj == Nx) && (kk == 1) %(N,N,1) , 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + + elseif (ii == 1) && (jj == Nx) && (kk == Ny) %(1,N,N), 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + elseif (ii == Nz) && (jj == 1) && (kk == Ny) %(N,1,N), 3 faces blocked + + n_check = 3; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk-1; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %Scenario when 1 side is blocked (moving along 1 face, in the middle/centre box, no edges), 6 cases + + + + elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == 1) %if on jj boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %jj boundary + + check_pos(1,1) = ii + 1; + check_pos(2,1) = ii - 1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj+1; %only option for jj to increase is to go up + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + + elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == 1) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %ii boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii+1; %only option for ii to increase is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == 1) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %kk boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + check_pos(5,1) = ii; %only option for ii to increase is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + check_pos(5,3) = kk+1; %only way is for kk to increase + + elseif (ii > 1) && (ii < Nz) && (kk >1) && (kk < Ny) && (jj == Nx) %if on jj boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %jj boundary + + check_pos(1,1) = ii + 1; + check_pos(2,1) = ii - 1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj-1; %only option for jj to decrease is to go up + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + + elseif (jj > 1) && (jj < Nx) && (kk > 1) && (kk < Ny) && (ii == Nz) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %ii boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + check_pos(5,1) = ii-1; %only option for ii to decrease is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + check_pos(5,3) = kk; + + elseif (jj > 1) && (jj < Nx) && (ii > 1) && (ii < Nz) && (kk == Ny) %if on ii boundary, one face blocked off + + n_check = 5; + check_pos = zeros(n_check, 5); + + %kk boundary + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + check_pos(5,1) = ii; %only option for ii to increase is to go up + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + check_pos(5,3) = kk-1; %only way is for kk to increase + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %Now: 2 sides blocked (edges of the faces), 12 cases + + %first 4: k is in middle + elseif (ii == 1) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + elseif (ii == Nz) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + elseif (ii == Nz) && (jj == 1) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii-1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + + elseif (ii == 1) && (jj == Nx) && (kk > 1) && (kk < Ny) %if on ii boundary, one face blocked off + + n_check = 4; + check_pos = zeros(n_check, 4); + + %kk boundary + + check_pos(1,1) = ii+1; + check_pos(2,1) = ii; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk+1; + check_pos(4,3) = kk-1; + + %next 4: i is in middle + + elseif (kk == 1) && (jj == 1) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (jj == Nx) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (jj == 1) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj+1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == 1) && (jj == Nx) && (ii > 1) && (ii < Nz) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + + check_pos(1,2) = jj; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + %next 4 cases: j is in middle + + elseif (kk == 1) && (ii == 1) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii+1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (ii == Nz) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii-1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == Ny) && (ii == 1) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii+1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk-1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + elseif (kk == 1) && (ii == Nz) && (jj > 1) && (jj < Nx) + + n_check = 4; + check_pos = zeros(n_check, 4); + + + check_pos(1,1) = ii; + check_pos(2,1) = ii-1; + check_pos(3,1) = ii; + check_pos(4,1) = ii; + + check_pos(1,2) = jj; + check_pos(2,2) = jj; + check_pos(3,2) = jj+1; + check_pos(4,2) = jj-1; + + check_pos(1,3) = kk+1; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %centre cubes - no sides are blocked, only need 1 case + + + elseif ((ii > 1) && (ii < Nz) && (jj >1) && (jj < Nx) && (kk > 1) && (kk < Ny)) + + + n_check = 6; + check_pos = zeros(n_check, 3); + + check_pos(1,1) = ii; + check_pos(2,1) = ii; + check_pos(3,1) = ii+1; + check_pos(4,1) = ii-1; + check_pos(5,1) = ii; + check_pos(6,1) = ii; + + check_pos(1,2) = jj+1; + check_pos(2,2) = jj-1; + check_pos(3,2) = jj; + check_pos(4,2) = jj; + check_pos(5,2) = jj; + check_pos(6,2) = jj; + + check_pos(1,3) = kk; + check_pos(2,3) = kk; + check_pos(3,3) = kk; + check_pos(4,3) = kk; + check_pos(5,3) = kk+1; + check_pos(6,3) = kk-1; + else + + disp("Error- cant find a situation") + + disp(ii) + disp(jj) + disp(kk) + + disp(Nz) + disp(Nx) + disp(Ny) + + end + + %check if values are outside of our domain + + %disp("checking for errors") + + for i = 1:n_check + if (check_pos(i,2) > Nx) + disp("i out of range") + + pause + end + + if (check_pos(i,1) > Nz) + disp("j out of range") + + pause + end + + + if (check_pos(i,3) > Ny) + disp("k out of range") + + pause + end + end + + + for i = 1:n_check + + %|| this means 'or' + % we are checking if we are at a boundary + + if (connected_mat(check_pos(i,1), check_pos(i,2), check_pos(i,3)) > 0) && (kr_phase_mat(ii,jj,kk) > tol) && (connected_mat(ii,jj,kk) == 0) + + connected_mat(ii,jj,kk) = 1; + n_changed = n_changed + 1; + end + + end + + + if (connected_mat(ii,jj,kk) == 1) && (ii == end_val) + connected = 1; + changing = 0; + break + else + connected = 0; + end + end + end + end + end + + if (n_changed == 0) + changing = 0; + end + + end +end +end diff --git a/Functions/check_box_connectivity.m b/src/check_box_connectivity.m similarity index 96% rename from Functions/check_box_connectivity.m rename to src/check_box_connectivity.m index c010ae8..522262e 100644 --- a/Functions/check_box_connectivity.m +++ b/src/check_box_connectivity.m @@ -1,72 +1,72 @@ -function [invaded_mat] = check_box_connectivity(pc, Pe_mat, Nz, Nx, Ny) - -invaded_mat = false(Nz, Nx, Ny); -changing = true; - -while (changing) - changing = false; - for ii = 1:Nz - if ii == 1 - near_check_z = 1; - elseif ii == Nz - near_check_z = -1; - else - near_check_z = [1,-1]; - end - - is_boundary_z = (ii == 1) || (ii == Nz); - - for jj = 1:Nx - - if jj == 1 - near_check_x = 1; - elseif jj == Nx - near_check_x = -1; - else - near_check_x = [1,-1]; - end - - is_boundary_x = (jj == 1) || (jj == Nx); - - for kk = 1:Ny - - if kk == 1 - near_check_y = 1; - elseif kk == Ny - near_check_y = -1; - else - near_check_y = [1,-1]; - end - - is_boundary_y = (kk == 1) || (kk == Ny); - - if pc <= Pe_mat(ii,jj,kk) - continue; - end - - if ~invaded_mat(ii,jj,kk) && (is_boundary_z || is_boundary_x || is_boundary_y) - invaded_mat(ii,jj,kk) = true; - changing = true; - end - - if invaded_mat(ii,jj,kk) - continue; - end - - for inc_z = near_check_z - for inc_x = near_check_x - for inc_y = near_check_y - if invaded_mat(ii+inc_z, jj+inc_x, kk+inc_y) - invaded_mat(ii,jj,kk) = true; - changing = true; - break; - end - end - end - end - - end - end - end - -end +function [invaded_mat] = check_box_connectivity(pc, Pe_mat, Nz, Nx, Ny) + +invaded_mat = false(Nz, Nx, Ny); +changing = true; + +while (changing) + changing = false; + for ii = 1:Nz + if ii == 1 + near_check_z = 1; + elseif ii == Nz + near_check_z = -1; + else + near_check_z = [1,-1]; + end + + is_boundary_z = (ii == 1) || (ii == Nz); + + for jj = 1:Nx + + if jj == 1 + near_check_x = 1; + elseif jj == Nx + near_check_x = -1; + else + near_check_x = [1,-1]; + end + + is_boundary_x = (jj == 1) || (jj == Nx); + + for kk = 1:Ny + + if kk == 1 + near_check_y = 1; + elseif kk == Ny + near_check_y = -1; + else + near_check_y = [1,-1]; + end + + is_boundary_y = (kk == 1) || (kk == Ny); + + if pc <= Pe_mat(ii,jj,kk) + continue; + end + + if ~invaded_mat(ii,jj,kk) && (is_boundary_z || is_boundary_x || is_boundary_y) + invaded_mat(ii,jj,kk) = true; + changing = true; + end + + if invaded_mat(ii,jj,kk) + continue; + end + + for inc_z = near_check_z + for inc_x = near_check_x + for inc_y = near_check_y + if invaded_mat(ii+inc_z, jj+inc_x, kk+inc_y) + invaded_mat(ii,jj,kk) = true; + changing = true; + break; + end + end + end + end + + end + end + end + +end diff --git a/src/downscale.m b/src/downscale.m new file mode 100644 index 0000000..6f3e388 --- /dev/null +++ b/src/downscale.m @@ -0,0 +1,53 @@ +function [subscale_dims, sub_porosity, sub_permeability, sub_entry_pressures, porosity] ... + = downscale(dr, params, options) + +subscale_dims = max(ceil(dr./options.subscale),3); +subscale_dims(mod(subscale_dims,2)==0) = subscale_dims(mod(subscale_dims,2)==0)+1; + +%% Permut arrays to satisfy (z,x,y) basis +idx_permut = [3,1,2]; +sub_dims_permut = subscale_dims(idx_permut); +dr_permut = dr(idx_permut); +corr_lens = params.poro.corr_lens(idx_permut); + +%% generate fine-scale porosity + +% TODO: consider using extreme value distribution `evrnd` +sub_porosity = rsgen3D(sub_dims_permut,corr_lens./dr_permut,@(dims) randn(dims).* params.poro.std_dev + params.poro.mean); +sub_porosity = max(sub_porosity,0); +porosity = sum(sub_porosity .* prod(dr_permut./double(sub_dims_permut)),'all')/prod(dr_permut); + +%% calculate fine-scale permeability + +sub_permeability = params.perm_corr(sub_porosity); + +%% calculate fine-scale sub_entry_pressures + +sub_entry_pressures = calc_entry_pressure(sub_porosity, sub_permeability, ... + params.capil.j_entry, params.capil.angle, params.capil.tension); + +end + +function [f] = rsgen3D(dims,corr_lens,dist) +arguments + dims (1,3) uint32 + corr_lens (1,3) double + dist = @(N) randn(N) +end + +z = linspace(-0.5,0.5,dims(1)); +x = linspace(-0.5,0.5,dims(2)); +y = linspace(-0.5,0.5,dims(3)); + +[X,Z,Y] = meshgrid(x,z,y); + +D = dist(dims.*2-1); + +% Gaussian filter +K = (X./corr_lens(2)).^2 + (Y./corr_lens(3)).^2 + (Z./corr_lens(1)).^2; +K = exp(K.*(-2)); +K = K./sum(K(:)); + +f = convn(D,K,'valid'); + +end diff --git a/src/downscale_upscale.m b/src/downscale_upscale.m new file mode 100644 index 0000000..b6aeda3 --- /dev/null +++ b/src/downscale_upscale.m @@ -0,0 +1,24 @@ +function [porosity, Kabs, sw_upscaled, pc_upscaled, krg, krw] = ... + downscale_upscale(dr, saturations, params, options) +arguments + dr (1,3) double + saturations (1,:) double + params (1,1) struct + options (1,1) struct +end + +% TODO: implement multi-scale approach for really small small_scale + +[subscale_dims, sub_porosity, sub_permeability, sub_entry_pressures, porosity] ... + = downscale(dr, params, options); + +if porosity <= 0 + return; +end + +[Kabs, sw_upscaled, pc_upscaled, krw, krg] = upscale(... + dr, saturations, params, options, ... + subscale_dims, sub_porosity, sub_permeability, sub_entry_pressures ... + ); + +end diff --git a/src/gen_options.m b/src/gen_options.m new file mode 100644 index 0000000..bb0e0a0 --- /dev/null +++ b/src/gen_options.m @@ -0,0 +1,9 @@ +function [options] = gen_options() +options = struct(... + 'subscale', 20*meter(), ... + 'sat_num_points', 16, ... + 'sat_min', 0.01, ... + 'sat_tol', 1e-8, ... + 'perm_min_threshold', 1e-8 * milli() * darcy() ... + ); +end diff --git a/src/gen_params.m b/src/gen_params.m new file mode 100644 index 0000000..635f76e --- /dev/null +++ b/src/gen_params.m @@ -0,0 +1,45 @@ +function [params] = gen_params() +params = struct(... + 'poro', struct(... + 'mean', 0.23,... + 'std_dev', 0.049,... + 'corr_lens',[200,20,0.1].*meter() ... + ),... + 'perm_corr', @(poro) 10.^(log10(poro).*4.19+5.17) * milli() * darcy(),... NOTE: produces milli Darcy + 'capil', struct(... + 'lambda', 1.75, ... + 'angle', deg2rad(37), ... + 'tension', 36 * milli() * newton() / meter(), ... + 'pres_entry', 10 * kilo() * pascal(), ... + 'sw_barrier', 0 ... + ),... + 'rel_perm', struct('sw_resid', 0.01)... + ); + +params.capil.pres_func = @(sw, entry_pressure) brooks_corey(sw, params.capil.sw_barrier, entry_pressure, params.capil.lambda); + +params.capil.pres_func_inv = @(pc, entry_pressure) brooks_corey_inv(pc, params.capil.sw_barrier, entry_pressure, params.capil.lambda); + +params.capil.pres_deriv = @(sw, entry_pressure) brooks_corey_deriv(sw, params.capil.sw_barrier, entry_pressure,params.capil.lambda); + +params.capil.j_entry = params.capil.pres_entry / (sqrt(params.poro.mean/params.perm_corr(params.poro.mean)) .* cos(params.capil.angle) .* params.capil.tension); + +params.rel_perm.wat_func = @(sw) max((sw - params.rel_perm.sw_resid)./(1 - params.rel_perm.sw_resid),0); +params.rel_perm.gas_func = @(sg) min(sg./(1 - params.rel_perm.sw_resid),1); +end + +function pc = brooks_corey(sw,sw_resid,p_entry,lambda) +if sw <= sw_resid + pc = Inf; + return; +end +sw = min(sw,1); +pc = p_entry .* ((1-sw_resid)./(sw-sw_resid)).^(1/lambda); +end + +function dpc_dsw = brooks_corey_deriv(sw,sw_resid,p_entry,lambda) +A = p_entry .* (1-sw_resid).^(1/lambda); +B = sw_resid; +C = -(1/lambda); +dpc_dsw = A.*C .* (sw - B).^(C-1); +end diff --git a/src/monotonize.m b/src/monotonize.m new file mode 100644 index 0000000..a01960f --- /dev/null +++ b/src/monotonize.m @@ -0,0 +1,8 @@ +function y = monotonize(x,y,direction) + for i = 2:length(x) + if (y(i) - y(i-1)) * direction >= 0 + continue; + end + y(i) = y(i-1); + end +end diff --git a/src/ogs/cell_to_curve.m b/src/ogs/cell_to_curve.m new file mode 100644 index 0000000..92695b6 --- /dev/null +++ b/src/ogs/cell_to_curve.m @@ -0,0 +1,7 @@ +function [curve_mapping, permeable_cells] = cell_to_curve(G,perm) + permeable_cells = 1:G.cells.num; + is_permeable = geomean(perm,2,"omitnan") > 0; + permeable_cells(is_permeable == 0) = []; + curve_mapping = zeros(prod(G.cartDims),1); + curve_mapping(G.cells.indexMap(permeable_cells)) = 1:length(permeable_cells); +end diff --git a/src/ogs/generate_sfn.m b/src/ogs/generate_sfn.m new file mode 100644 index 0000000..50cabc4 --- /dev/null +++ b/src/ogs/generate_sfn.m @@ -0,0 +1,73 @@ +function [file_name] = generate_sfn(mask, saturations, pres_upscaled, krw, krg,prefix,inc_ext) +arguments + mask (:,1) logical + saturations (1,:) double + pres_upscaled (:,:) double + krw (:,3,:) double + krg (:,3,:) double + prefix char + inc_ext char = '.inc' +end +file_name = [prefix,'sfn',inc_ext]; + +file_id = fopen(file_name,'w','native','UTF-8'); + + +for direction=1:3 + write_tables_for_direction(mask,file_id,krw, krg,saturations,pres_upscaled,direction); +end +fclose(file_id); +end + +function write_tables_for_direction(mask,file_id, krw, krg,saturations,pres_upscaled,direction) +for cell_index = 1:length(mask) + if ~mask(cell_index) + continue; + end + + table_num = cell_index + length(mask) * (direction-1); + + fprintf(file_id,'CHARACTERISTIC_CURVES sfn_%u\n',table_num); + fprintf(file_id,'%s\n%s\n',... + 'TABLE swfn_table',... + 'PRESSURE_UNITS Bar'); + + sw = saturations; + write_sfn_table(file_id,"SWFN",... + sw, ... + squeeze(krw(cell_index,direction,:)), ... + pres_upscaled(cell_index,:)); + + fprintf(file_id,'%s\n%s\n%s\n',... + 'END',... + 'TABLE sgfn_table',... + 'PRESSURE_UNITS Bar'); + + sg = 1 - sw; + sg = sg(end:-1:1); + write_sfn_table(file_id,"SGFN",... + sg, ... + squeeze(krg(cell_index,direction,end:-1:1)), ... + zeros(size(sg))); + + fprintf(file_id,'%s\n%s\n',... + 'END',... + 'END'); +end +end + +function write_sfn_table(file_id,table_name,sat,kr,pc) +arguments + file_id + table_name string + sat (1,:) double + kr (1,:) double + pc (1,:) double +end + +data = [sat;kr;pc]; + +fprintf(file_id,'%s\n',table_name); +fprintf(file_id,'%f\t%f\t%f\n',data); +fprintf(file_id,'%s\n','/'); +end diff --git a/src/ogs/write_curve_nums.m b/src/ogs/write_curve_nums.m new file mode 100644 index 0000000..391dd07 --- /dev/null +++ b/src/ogs/write_curve_nums.m @@ -0,0 +1,12 @@ +function write_curve_nums(file_name,keyword,curve_mapping, fallback_curve_num, offset) +line_ending = "\n"; + +writelines(keyword,file_name,LineEnding=line_ending,WriteMode='overwrite'); + +mapping_to_write = curve_mapping; +mapping_to_write(curve_mapping==0) = fallback_curve_num; +mapping_to_write(curve_mapping~=0) = mapping_to_write(curve_mapping~=0) + offset; +writematrix(mapping_to_write,file_name,FileType="text",WriteMode="append",Delimiter='space'); + +writelines(["/",""],file_name,LineEnding=line_ending,WriteMode='append'); +end diff --git a/src/ogs/write_mappings.m b/src/ogs/write_mappings.m new file mode 100644 index 0000000..c0180f3 --- /dev/null +++ b/src/ogs/write_mappings.m @@ -0,0 +1,30 @@ +function write_mappings(prefix,G,perm,poro) +[curve_mapping, permeable_cells] = cell_to_curve(G, perm); + +keywords = ["KRNUMX","KRNUMY","KRNUMZ","SATNUM"]; + +num_permeable_cells = length(permeable_cells); +num_offset = [0,num_permeable_cells,num_permeable_cells*2,0]; + +parfor keyword_num = 1:length(keywords) + keyword = keywords(keyword_num); + file_name = join([prefix,keyword,".inc"],''); + write_curve_nums(file_name,keyword,curve_mapping,1,num_offset(keyword_num)); +end + +perm_mapping = zeros(G.cells.num,3); +perm_mapping(G.cells.indexMap,:) = perm; +keywords = ["PERMX","PERMY","PERMZ"]; +parfor keyword_num = 1:length(keywords) + keyword = keywords(keyword_num); + file_name = join([prefix,keyword,".grdecl"],''); + write_curve_nums(file_name,keyword,perm_mapping(:,keyword_num),0,0); +end + +poro_mapping = zeros(G.cells.num,1); +poro_mapping(G.cells.indexMap) = poro; +keyword = "PORO"; +file_name = join([prefix,keyword,".grdecl"],''); +write_curve_nums(file_name,"PORO",poro_mapping,0,0); + +end diff --git a/src/units/barsa.m b/src/units/barsa.m new file mode 100644 index 0000000..92e0aef --- /dev/null +++ b/src/units/barsa.m @@ -0,0 +1,4 @@ +function [mult] = barsa() +mult = 1e5; +end + diff --git a/src/units/darcy.m b/src/units/darcy.m new file mode 100644 index 0000000..6a3c76d --- /dev/null +++ b/src/units/darcy.m @@ -0,0 +1,4 @@ +function [mult] = darcy() +mult = 9.869232667160130e-13 * meter()^2; +end + diff --git a/src/units/kilo.m b/src/units/kilo.m new file mode 100644 index 0000000..0e807ca --- /dev/null +++ b/src/units/kilo.m @@ -0,0 +1,4 @@ +function [mult] = kilo() +mult = 1e3; +end + diff --git a/src/units/meter.m b/src/units/meter.m new file mode 100644 index 0000000..bcf4f0f --- /dev/null +++ b/src/units/meter.m @@ -0,0 +1,4 @@ +function [mult] = meter() +mult = 1.0; +end + diff --git a/src/units/milli.m b/src/units/milli.m new file mode 100644 index 0000000..97ac8f5 --- /dev/null +++ b/src/units/milli.m @@ -0,0 +1,4 @@ +function [mult] = milli() +mult = 1e-3; +end + diff --git a/src/units/newton.m b/src/units/newton.m new file mode 100644 index 0000000..ef5e4b8 --- /dev/null +++ b/src/units/newton.m @@ -0,0 +1,3 @@ +function [mult] = newton() +mult = 1.0; +end \ No newline at end of file diff --git a/src/units/pascal.m b/src/units/pascal.m new file mode 100644 index 0000000..f6f613c --- /dev/null +++ b/src/units/pascal.m @@ -0,0 +1,4 @@ +function [mult] = pascal() +mult = 1.0; +end + diff --git a/src/upscale.m b/src/upscale.m new file mode 100644 index 0000000..293d1ff --- /dev/null +++ b/src/upscale.m @@ -0,0 +1,292 @@ +function [Kabs, sw_upscaled, pc_upscaled, krw, krg] = upscale(... + dr, saturations, params, options, ... + downscale_dims, porosities, permeabilities, entry_pressures) + +if max(porosities,[],'all') <= 0 + error('inactive cell'); +end + +dr_sub = dr ./ downscale_dims; + +Kabs = Calculate_Kabs(permeabilities, dr_sub(1),dr_sub(2),dr_sub(3), 1, 1, 1); + +sw_upscaled = saturations; +pc_upscaled = zeros(size(saturations)); + +Nz_sub = downscale_dims(3); +Nx_sub = downscale_dims(1); +Ny_sub = downscale_dims(2); + +krg = zeros(3,length(sw_upscaled)); +krw = zeros(3,length(sw_upscaled)); + +tol_kr = options.perm_min_threshold; + +for index_saturation = 1:length(saturations) + + sw_target = saturations(index_saturation); + + pc_mid = params.capil.pres_func(sw_target, mean(entry_pressures,'all')); + sw_mid = sw_target; + for iteration_num=1:1000 + + [pc_mid_tot, sw_mid, pc_mid, invaded_mat_mid, converged] = mip_iteration(... + sw_target, dr, entry_pressures, porosities, pc_mid, ... + Nz_sub, Nx_sub, Ny_sub,... + params.capil.pres_func_inv, params.capil.pres_deriv,... + options.sat_tol ... + ); + + if converged + break; + end + end + + sw_upscaled(index_saturation) = sw_mid; + pc_upscaled(index_saturation) = pc_mid_tot; + + sw = invaded_mat_mid .* sw_mid + ~invaded_mat_mid .* 1; + + kg_mat_local = params.rel_perm.gas_func(1-sw); + kw_mat_local = params.rel_perm.wat_func(sw); + + kg_mat_local = kg_mat_local.*permeabilities; % FIXME scaling to absolute permeabilities + kw_mat_local = kw_mat_local.*permeabilities; % FIXME scaling to absolute permeabilities + + kg_phase_connected_local = zeros(1,3); + kw_phase_connected_local = zeros(1,3); + + [kg_phase_connected_local(1), ~] = check_axis_connectivity(kg_mat_local, Nx_sub, Nz_sub,Ny_sub,tol_kr, 1, Ny_sub); %z + [kg_phase_connected_local(2), ~] = check_axis_connectivity(kg_mat_local, Nx_sub, Nz_sub,Ny_sub,tol_kr, 2, Nx_sub); %y + [kg_phase_connected_local(3), ~] = check_axis_connectivity(kg_mat_local, Nx_sub, Nz_sub,Ny_sub,tol_kr, 3, Nz_sub); %x + + [kw_phase_connected_local(1), ~] = check_axis_connectivity(kw_mat_local, Nx_sub, Nz_sub,Ny_sub,tol_kr, 1, Ny_sub);%z + [kw_phase_connected_local(2), ~] = check_axis_connectivity(kw_mat_local, Nx_sub, Nz_sub,Ny_sub,tol_kr, 2, Nx_sub);%y + [kw_phase_connected_local(3), ~] = check_axis_connectivity(kw_mat_local, Nx_sub, Nz_sub,Ny_sub,tol_kr, 3, Nz_sub);%x + + + [krg(:,index_saturation), krw(:,index_saturation)] = calc_phase_permeabilities(... + dr_sub, Kabs, ... + kg_phase_connected_local, kw_phase_connected_local,... + kg_mat_local, kw_mat_local... + ); +end + +sw_upscaled(end) = 1; +[sw_upscaled,unique_idx] = unique(sw_upscaled); +pc_upscaled = pc_upscaled(unique_idx); +krg = krg(:,unique_idx); +krw = krw(:,unique_idx); + +[sw_upscaled,sorted_idx] = sort(sw_upscaled); +pc_upscaled = pc_upscaled(sorted_idx); +krg = krg(:,sorted_idx); +krw = krw(:,sorted_idx); + +pc_upscaled(isinf(pc_upscaled)) = max(pc_upscaled(~isinf(pc_upscaled))); + +if ~allfinite(pc_upscaled) + disp(pc_upscaled); +end + +end + +function [pc_mid_tot, sw_mid, pc_mid, invaded_mat_mid, converged] = mip_iteration(... + sw_target, dr, entry_pressures, porosities, pc_mid,... + Nz_sub, Nx_sub, Ny_sub, ... + pc_func_inv, pc_deriv, ... + tol_sw) + +invaded_mat_mid = check_box_connectivity(pc_mid, entry_pressures, Nz_sub, Nx_sub, Ny_sub); + +volume = prod(dr); +sub_volume = volume./double(Nz_sub*Nx_sub*Ny_sub); +pore_volumes = porosities .* sub_volume; +pore_volume = sum(pore_volumes,'all'); + +sub_sw_mid = invaded_mat_mid .* pc_func_inv(pc_mid,entry_pressures) + ~invaded_mat_mid .* 1; +sub_sw_mid(~isfinite(sub_sw_mid)) = 1; +sw_mid = sum(sub_sw_mid.*pore_volumes,'all')/pore_volume; + +pc_mid_tot = sum((1-sub_sw_mid).*pore_volumes.*pc_mid,"all")/(pore_volume*(1-sw_mid)); + +if sw_mid >=1 + pc_mid_tot = pc_mid; % FIXME remove +end + +sw_err = sw_target - sw_mid; +err = abs(sw_err); +converged = err <= tol_sw; +if converged + return; +end + +deriv = pc_deriv(sw_mid, mean(entry_pressures,'all')); + +dpc = sw_err*deriv; + +pc_mid = pc_mid + dpc * 0.5; +if ~isfinite(pc_mid) + error('') +end +pc_mid = max(pc_mid,min(entry_pressures,[],'all')); +end + +function [krg, krw] = calc_phase_permeabilities(... + dr_sub, Kabs,... + kg_phase_connected, kw_phase_connected, kg_mat, kw_mat ... + ) + +Kx = Kabs(1); +Ky = Kabs(2); +Kz = Kabs(3); + +Kalli = zeros(1,6); + +h1g = (kg_phase_connected(3) > 0); +h1w = (kw_phase_connected(3) > 0); +h2g = (kg_phase_connected(2) > 0); +h2w = (kw_phase_connected(2) > 0); +v1g = (kg_phase_connected(1) > 0); +v1w = (kw_phase_connected(1) > 0); +if h1g || h2g || v1g + Estimate = Calculate_Kabs(kg_mat, dr_sub(1), dr_sub(2), dr_sub(3), h1g, h2g, v1g); + Kalli(1) = Estimate(1); + Kalli(2) = Estimate(2); + Kalli(3) = Estimate(3); +end +if h1w || h2w || v1w + Estimate = Calculate_Kabs(kw_mat,dr_sub(1), dr_sub(2), dr_sub(3), h1w, h2w, v1w); + Kalli(4) = Estimate(1); + Kalli(5) = Estimate(2); + Kalli(6) = Estimate(3); +end + +Kalli([1,4]) = Kalli([1,4]) ./ Kx; +Kalli([2,5]) = Kalli([2,5]) ./ Ky; +Kalli([3,6]) = Kalli([3,6]) ./ Kz; + +krg = Kalli(1:3)'; +krw = Kalli(4:6)'; + +end + +function [Estimate] = Calculate_Kabs(Kai0,Lx,Ly,Lz,h1,h2,v1) % TODO support anisotropic inputs +Kh1 = 0; +Kh2 = 0; +Kv1 = 0; + +% effective simulation domain, along y direction +[x0,y0,z0] = size(Kai0); +Q = 1 * 10^-6 / 60; % (ml/min)->m3/s % FIXME: replace with input parameter +Muw = 1 * 10^-3; % cp-PaS % FIXME: replace with input parameter +Pout = 1.0; % FIXME: replace with input parameter + +x=x0+2; %for layer +y=y0+2; %for layer +z=z0+2; + +Kai = zeros(x,y,z); +Kai(2:end-1,2:end-1,2:end-1) = Kai0; +%no value for first and endlayer + +Condx = Kai(:) .* (Ly*Lz) ./ Muw ./Lx; +Condy = Kai(:) .* (Lx*Lz) ./ Muw ./Ly; +Condz = Kai(:) .* (Lx*Ly) ./ Muw ./Lz; + +Loc=zeros(x,y,z); +Ind = sub2ind(size(Loc),find(Condz>0)); +Loc(Ind) = (1:length(Ind)); +n = length(Ind) + 1; %+1 for constant Q + +Dim = find(Loc>0); +Dim1 = Dim - 1; +Dim2 = Dim + 1; +Dim3 = Dim - x; +Dim4 = Dim + x; +Dim5 = Dim - x*y; +Dim6 = Dim + x*y; +%Loc(Dim) already > 0 +A0 = zeros(length(Dim)+1,3); +A0(:,1) = [Loc(Dim);n]; +A0(:,2) = [Loc(Dim);n]; +ind = find(Loc(Dim1)>0); +A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2./(1./Condx(Dim(ind))+1./Condx(Dim1(ind))); +A1 = [Loc(Dim(ind)),Loc(Dim1(ind)), -2./(1./Condx(Dim(ind))+1./Condx(Dim1(ind)))]; +ind = find(Loc(Dim2)>0); +A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2./(1./Condx(Dim(ind))+1./Condx(Dim2(ind))); +A2 = [Loc(Dim(ind)),Loc(Dim2(ind)), -2./(1./Condx(Dim(ind))+1./Condx(Dim2(ind)))]; +ind = find(Loc(Dim3)>0); +A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2./(1./Condy(Dim(ind))+1./Condy(Dim3(ind))); +A3 = [Loc(Dim(ind)),Loc(Dim3(ind)), -2./(1./Condy(Dim(ind))+1./Condy(Dim3(ind)))]; +ind = find(Loc(Dim4)>0); +A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2./(1./Condy(Dim(ind))+1./Condy(Dim4(ind))); +A4 = [Loc(Dim(ind)),Loc(Dim4(ind)), -2./(1./Condy(Dim(ind))+1./Condy(Dim4(ind)))]; +ind = find(Loc(Dim5)>0); +A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2./(1./Condz(Dim(ind))+1./Condz(Dim5(ind))); +A5 = [Loc(Dim(ind)),Loc(Dim5(ind)), -2./(1./Condz(Dim(ind))+1./Condz(Dim5(ind)))]; +ind = find(Loc(Dim6)>0); +A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2./(1./Condz(Dim(ind))+1./Condz(Dim6(ind))); +A6 = [Loc(Dim(ind)),Loc(Dim6(ind)), -2./(1./Condz(Dim(ind))+1./Condz(Dim6(ind)))]; +A0b = A0; +%boundary y direction +if h1 + ind = find(Loc(Dim3)<=0); + A0(n,3) = A0(n,3) + 2.0.*sum(Condy(Dim(ind))); + A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2.0.*Condy(Dim(ind)); + Ain = [[Loc(Dim(ind));ind*0+n],[ind*0+n;Loc(Dim(ind))],[-2.0.*Condy(Dim(ind));-2.0.*Condy(Dim(ind))]]; + ind = find(Loc(Dim4)<=0); + A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2.0.*Condy(Dim(ind)); + Aall=[A0;A1;A2;A3;A4;A5;A6;Ain]; + A = sparse(Aall(:,1),Aall(:,2),Aall(:,3),n,n); + B = zeros(n,1); + B(Loc(Dim(ind))) = Pout * 2.0.*Condy(Dim(ind)); + B(end) = Q * 10^12; %due to D + X1=mldivide(A,B); + Pin = X1(end); + Ae = x0*z0*Lx*Lz; %m2 + Le = (y0+0)*Ly; %m,+0 considering bouyndary is half, +1 for one + Kh1 = Q/((Pin-Pout)*Ae/Le/Muw) * 10^12; %Darcy +end +if h2 + A0 = A0b; + ind = find(Loc(Dim5)<=0); + A0(n,3) = A0(n,3) + 2.0.*sum(Condz(Dim(ind))); + A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2.0.*Condz(Dim(ind)); + Ain = [[Loc(Dim(ind));ind*0+n],[ind*0+n;Loc(Dim(ind))],[-2.0.*Condz(Dim(ind));-2.0.*Condz(Dim(ind))]]; + ind = find(Loc(Dim6)<=0); + A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2.0.*Condz(Dim(ind)); + Aall=[A0;A1;A2;A3;A4;A5;A6;Ain]; + A = sparse(Aall(:,1),Aall(:,2),Aall(:,3),n,n); + B = zeros(n,1); + B(Loc(Dim(ind))) = Pout * 2.0.*Condz(Dim(ind)); + B(end) = Q * 10^12; %due to D + X1=mldivide(A,B); + Pin = X1(end); + Ae = x0*y0*Lx*Ly; %m2 + Le = (z0+0)*Lz; %m + Kh2 = Q/((Pin-Pout)*Ae/Le/Muw) * 10^12; %Darcy +end +if v1 + A0 = A0b; + ind = find(Loc(Dim1)<=0); + A0(n,3) = A0(n,3) + 2.0.*sum(Condx(Dim(ind))); + A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2.0.*Condx(Dim(ind)); + Ain = [[Loc(Dim(ind));ind*0+n],[ind*0+n;Loc(Dim(ind))],[-2.0.*Condx(Dim(ind));-2.0.*Condx(Dim(ind))]]; + ind = find(Loc(Dim2)<=0); + A0(Loc(Dim(ind)),3) = A0(Loc(Dim(ind)),3) + 2.0.*Condx(Dim(ind)); + Aall=[A0;A1;A2;A3;A4;A5;A6;Ain]; + A = sparse(Aall(:,1),Aall(:,2),Aall(:,3),n,n); + B = zeros(n,1); + B(Loc(Dim(ind))) = Pout * 2.0.*Condx(Dim(ind)); + B(end) = Q * 10^12; %due to D + X1=mldivide(A,B); + Pin = X1(end); + Ae = y0*z0*Ly*Lz; %m2 + Le = (x0+0)*Lx; %m + Kv1 = Q/((Pin-Pout)*Ae/Le/Muw) * 10^12; %Darcy +end + +Estimate = [Kh1,Kh2,Kv1]; + +end diff --git a/startup.m b/startup.m new file mode 100644 index 0000000..7ecff93 --- /dev/null +++ b/startup.m @@ -0,0 +1 @@ +addpath(genpath("./src")); diff --git a/strata_trapper.m b/strata_trapper.m new file mode 100644 index 0000000..ac789de --- /dev/null +++ b/strata_trapper.m @@ -0,0 +1,120 @@ +function strata_trapped = strata_trapper(G, mask, params, options, enable_waitbar) +arguments + G (1,1) struct + mask (:,1) logical + params (1,1) struct + options (1,1) struct + enable_waitbar (1,1) logical = false; +end + +perm_upscaled = zeros(G.cells.num, 3); +poro_upscaled = zeros(G.cells.num, 1); + +saturations = unique([linspace(options.sat_min,1,options.sat_num_points),params.rel_perm.sw_resid]); + +cap_pres_upscaled = zeros(G.cells.num,length(saturations)); +krw = zeros(G.cells.num, 3,length(saturations)); +krg = zeros(G.cells.num,3,length(saturations)); + +cells_num = min(length(mask),G.cells.num); + +if enable_waitbar + wb_queue = parallel.pool.DataQueue; + parforWaitbar(0,sum(mask)); + afterEach(wb_queue,@parforWaitbar); +end + +DR = [G.DX,G.DY,G.DZ]; + +parfor cell_index = 1:cells_num + if ~mask(cell_index) + continue; + end + + [porosity, Kabs, sw_upscaled, pc_upscaled, krg_cell, krw_cell] = ... + downscale_upscale(DR(cell_index,:), saturations , params, options); + + if porosity <= 0 + error('Resulted in non-positive upscaled porosity') + end + + poro_upscaled(cell_index) = porosity; + perm_upscaled(cell_index,:) = Kabs; + cap_pres_upscaled(cell_index,:) = interp1(sw_upscaled,pc_upscaled,saturations,"linear","extrap"); + + for i = 1:3 + krw_cell(i,:) = monotonize(sw_upscaled, krw_cell(i,:), 1); + krg_cell(i,:) = monotonize(sw_upscaled, krg_cell(i,:), -1); + end + + krw(cell_index,:,:) = interp1(sw_upscaled, krw_cell', saturations, "linear")'; + krg(cell_index,:,:) = interp1(sw_upscaled, krg_cell', saturations, "linear")'; + + if enable_waitbar + send(wb_queue,cell_index); + end +end + +krw(:,:,saturations<=params.rel_perm.sw_resid) = 0; +krg(:,:,saturations<=params.rel_perm.sw_resid) = 1; +krg(:,:,end) = 0; + +strata_trapped = struct(... + 'porosity', poro_upscaled, ... + 'permeability', perm_upscaled, ... + 'saturation', saturations,... + 'capillary_pressure', cap_pres_upscaled, ... + 'rel_perm_wat', krw, ... + 'rel_perm_gas', krg ... + ); + +if enable_waitbar + parforWaitbar(0,0,'ready'); +end + +end + +function parforWaitbar(~,max_iterations,~) +persistent state wb final_state timer_val last_reported_state last_reported_time + +if nargin == 2 + state = 0; + final_state = max_iterations; + wb = waitbar(state,sprintf('0/%u iterations', final_state),'Name','StrataTrapper'); + timer_val = tic(); + + last_reported_state = state; + last_reported_time = timer_val; + return; +end + +if ~isvalid(wb) + return; +end + +if nargin == 3 + elapsed = duration(seconds(toc(timer_val)),'Format','hh:mm:ss'); + message = sprintf('%u coarse cells in %s',final_state,elapsed); + waitbar(1,wb,message); + return; +end + +state = state + 1; +elapsed = toc(timer_val); + +time_to_report = (elapsed - last_reported_time) > 1; +state_to_report = (state - last_reported_state) > final_state * 0.001; + +if and(~time_to_report,~state_to_report) + return; +end + +last_reported_state = state; +last_reported_time = elapsed; + +t_expected = elapsed * final_state / state; +eta = duration(seconds(t_expected - elapsed),'Format','hh:mm:ss'); +elapsed = duration(seconds(elapsed),'Format','hh:mm:ss'); +message = sprintf('%u/%u coarse cells\n passed: %s | ETA: %s',state,final_state,elapsed,eta); +waitbar(state/final_state,wb,message); +end