diff --git a/@Ternary_model/create_mesh.m b/@Ternary_model/create_mesh.m index 2cba8028bfb160153bfc0aa6f4bebabe69f3dd8e..3fa8a6a811a598ae3ccd50bf23a521d899247fd8 100644 --- a/@Ternary_model/create_mesh.m +++ b/@Ternary_model/create_mesh.m @@ -2,36 +2,32 @@ function create_mesh(T) %CREATE_MESH creates mesh for solving Ternary_model. % Mesh density larger for steeper gradients. -if strcmp(T.mode, 'Constituent') % Start with left side of mesh, at r=0/x=0, with a given step size - x = linspace(0, 0.75, 40); - x = [x, x(end)+0.001/T.precision]; - while x(end)<-T.a - x_temp = x(end)-x(end-1); - phi_nm1 = T.phi_tot(x(end), T.a, T.b, T.e, T.u0, 0); - phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0); - while phi_nm1-phi_n > 0.002/T.precision +x = [0, 0.01/T.precision];%/10 +x = [0, 0.01/T.precision/10];%/10 + +while x<-2*T.a + x_temp = x(end)-x(end-1); + phi_nm1 = T.phi_tot(x(end), T.a, T.b, T.e, T.u0, 0); + phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0); + if x(end) <= -T.a + while phi_nm1-phi_n > 0.0002/T.precision x_temp = x_temp/2; phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0); end - x = [x, x(end)+x_temp]; - end - x = [x, x(end)+cumsum(fliplr(diff(x)))]; - x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 100*T.precision); - x_right = linspace(-10*T.a+0.1, T.system_size, 100*T.precision); - T.x = [x, x_middle, x_right]; - T.x = T.x(1:find(T.x-T.system_size>=0, 1)-1); - - % Fine mesh close to x0 if ic = 'Gauss' - [~, ind] = min(abs(T.x-T.x0)); - if (T.x(ind+20)-T.x(ind) > 0.1) && strcmp(T.ic, 'Gauss') - T.x = [T.x(1:ind-21), linspace(T.x(ind-20), T.x(ind+20),... - ceil((T.x(ind+20)-T.x(ind-20))*200)),... - T.x(ind+21:end)]; + elseif x(end) > -T.a + while (phi_nm1-phi_n < 0.0001/T.precision) && ... + (x_temp < 0.01/T.precision) + x_temp = 2*x_temp; + phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0); + end end -% T.x = [linspace(0, 0.99, 100), linspace(0.9901, 1.01, 2000), linspace(1.0161, 2, 100),... -% linspace(2.01, 10, 100), linspace(10.01, 300, 200)]; -elseif strcmp(T.mode, 'Client') - T.x = linspace(0, 20, 24000); + x = [x, x(end)+x_temp]; end +x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 200*T.precision); +x_right = linspace(-10*T.a+0.1, T.system_size, 200*T.precision); +x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 600*T.precision); +x_right = linspace(-10*T.a+0.1, T.system_size, 600*T.precision); +T.x = [x, x_middle, x_right]; +T.x = T.x(1:find(T.x-T.system_size>=0, 1)-1); end