diff options
Diffstat (limited to 'CSXCAD/matlab/SmoothMeshLines2.m')
| -rw-r--r-- | CSXCAD/matlab/SmoothMeshLines2.m | 318 |
1 files changed, 318 insertions, 0 deletions
diff --git a/CSXCAD/matlab/SmoothMeshLines2.m b/CSXCAD/matlab/SmoothMeshLines2.m new file mode 100644 index 0000000..5d347c0 --- /dev/null +++ b/CSXCAD/matlab/SmoothMeshLines2.m @@ -0,0 +1,318 @@ +function lines = SmoothMeshLines2(lines, max_res, ratio, varargin) +%lines = SmoothMeshLines2( lines, max_res [, ratio, varargin] ) +% +% Create smooth mesh lines. +% +% input: +% lines: 1xn vector of (fixed) mesh lines +% max_res: maximum distance between any two lines (e.g. lambda/10) +% ratio: (optional) default: 1.3 +% +% optional variable arguments ('key', value) +% CheckMesh: Do a final mesh check (default is true) +% allowed_max_ratio: allow only a given max. grading ratio +% (default --> ratio*1.25) +% output: +% lines: 1xn vector of (smoothed) mesh lines +% +% example: +% mesh.x = [0 100 2300 2400]; +% mesh.x = SmoothMeshLines2( mesh.x, 43 ); +% +% todo: +% - if gaps cannot be completely filled because of the ratio restriction, +% the spacing is not optimal. SmoothRange() needs to be optimized. +% - SmoothRange() has special handling if stop_taper is to long; can this +% also happen for start_taper? +% +% CSXCAD matlab interface +% ----------------------- +% (C) 2010 Sebastian Held <sebastian.held@uni-due.de> +% See also SmoothMeshLines + +if (numel(lines)<2) + return +end + +if (nargin<3) + ratio = 1.3; +end + +check_mesh = true; +max_ratio = ratio*1.25; + +for vn=1:2:numel(varargin) + if (strcmpi(varargin{vn},'CheckMesh')) + check_mesh = varargin{vn+1}; + end + if (strcmpi(varargin{vn},'allowed_max_ratio')) + max_ratio = varargin{vn+1}; + end +end + +% +% | | | | +% | gap(1) | gap(2) | ... | +% | | | | +% +% lines(1) lines(2) ... lines(end) + +lines = unique(sort(lines)); + +for n=1:(numel(lines)-1) + old_gap(n).start_res = -1; + old_gap(n).stop_res = -1; + old_gap(n).lines = []; +end + +if numel(lines) == 2 + % special case + addLines = SmoothRange( lines(1), lines(2), [], [], max_res, ratio ); + lines = sort(unique([lines addLines])); + [EC pos] = CheckMesh(lines,0,max_res,ratio); + return +end + +while 1 + gap = calc_gaps(lines,old_gap,max_res); + + % determine gap to process + index = find( [old_gap(2:end).start_res] ~= [gap(2:end).start_res], 1, 'first' ) + 1; + if isempty(index) + index = find( [old_gap(1:end-1).stop_res] ~= [gap(1:end-1).stop_res], 1, 'first' ); + end + if isempty(index) + break; % done + end + + start_res = min( max_res, calc_start_res(index,lines,old_gap) ); + stop_res = min( max_res, calc_stop_res(index,lines,old_gap) ); + + % create new lines + old_gap(index).lines = SmoothRange( lines(index), lines(index+1), start_res, stop_res, max_res, ratio ); + + % remember gap setting + old_gap(index).start_res = start_res; + old_gap(index).stop_res = stop_res; + + % debug +% plot_lines +end + +% merge lines +for n=1:numel(old_gap) + lines = [lines old_gap(n).lines]; +end +lines = sort(unique(lines)); + +if (check_mesh) + CheckMesh(lines,0,max_res,max_ratio,0); +end + +end % SmoothMeshLines2() + +%% ------------------------------------------------------------------------ +function gap = calc_gaps(lines,old_gap,max_res) + temp_lines = lines; + for n=1:numel(old_gap) + temp_lines = [temp_lines old_gap(n).lines]; + end + temp_lines = sort(unique(temp_lines)); + + if numel(temp_lines) == 2 + gap(1).start_res = inf; % resolution not fixed + gap(1).stop_res = max_res; + return + end + + gap(1).start_res = inf; % resolution not fixed + idx = interp1( temp_lines, 1:numel(temp_lines), lines(2), 'nearest' ); + gap(1).stop_res = min( max_res, temp_lines(idx+1) - temp_lines(idx) ); + for n=2:numel(lines)-2 + idx = interp1( temp_lines, 1:numel(temp_lines), lines(n), 'nearest' ); + gap(n).start_res = min( max_res, temp_lines(idx) - temp_lines(idx-1) ); + idx = interp1( temp_lines, 1:numel(temp_lines), lines(n+1), 'nearest' ); + gap(n).stop_res = min( max_res, temp_lines(idx+1) - temp_lines(idx) ); + gap(n).lines = old_gap(n).lines; + end + idx = interp1( temp_lines, 1:numel(temp_lines), lines(end-1), 'nearest' ); + gap(numel(lines)-1).start_res = min( max_res, temp_lines(idx) - temp_lines(idx-1) ); + gap(numel(lines)-1).stop_res = inf; % resolution not fixed +end % calc_gaps() + + +%% ------------------------------------------------------------------------ +function plot_lines + temp_lines = lines; + for n=1:numel(old_gap) + temp_lines = [temp_lines old_gap(n).lines]; + end + temp_lines = sort(unique(temp_lines)); + + plot( temp_lines, ones(size(temp_lines)), 'r+' ); + hold on + plot( lines, ones(size(lines)), 'bo' ); + hold off +end % plot_lines + + +%% ------------------------------------------------------------------------ +function res = calc_start_res(pos,lines,old_gap) + if (pos < 2) || (pos > numel(lines)) + res = []; + return + end + temp_lines = lines; + for n=1:numel(old_gap) + temp_lines = [temp_lines old_gap(n).lines]; + end + temp_lines = sort(unique(temp_lines)); + + idx = interp1( temp_lines, 1:numel(temp_lines), lines(pos), 'nearest' ); + res = temp_lines(idx) - temp_lines(idx-1); +end % calc_res() + + +%% ------------------------------------------------------------------------ +function res = calc_stop_res(pos,lines,old_gap) + if (pos < 1) || (pos >= numel(lines)-1) + res = []; + return + end + temp_lines = lines; + for n=1:numel(old_gap) + temp_lines = [temp_lines old_gap(n).lines]; + end + temp_lines = sort(unique(temp_lines)); + + idx = interp1( temp_lines, 1:numel(temp_lines), lines(pos+1), 'nearest' ); + res = temp_lines(idx+1) - temp_lines(idx); +end % calc_res() + + +%% ------------------------------------------------------------------------ +function lines = SmoothRange(start, stop, start_res, stop_res, max_res, ratio) + +if (nargin<6) + ratio = 1.3; +end + +if isempty(start_res) && isempty(stop_res) + % special case: fill entire range with max_res + n1 = ceil( (stop-start) / max_res ) + 1; + lines = linspace( start, stop, n1 ); + return +end + +if isempty(start_res) + % special case: taper from stop_res at stop to max_res at start (if + % possible) + taper = stop_res*ratio; + stop_taper = stop; + while (taper*ratio < max_res) && (stop_taper(1)-taper > start) + stop_taper = [stop_taper(1)-taper stop_taper]; + taper = taper*ratio; + end + if (stop_taper(1) - start) >= max_res + % fill with equidistant lines + n1 = ceil( (stop_taper(1)-start) / max_res ) + 1; + stop_taper = [linspace(start,stop_taper(1),n1) stop_taper]; + else + % not enough space for entire taper + stop_taper(1) = []; % likely too near to start + if numel(stop_taper) > 0 + stop_taper = [(start+stop_taper(1))/2 stop_taper]; % create a centered line + end + end + lines = sort(unique(stop_taper(stop_taper>=start))); + return +end + +if isempty(stop_res) + % special case: taper from stop_res at stop to max_res at start (if + % possible) + taper = start_res*ratio; + start_taper = start; + while (taper*ratio < max_res) && (start_taper(end)+taper < stop) + start_taper = [start_taper start_taper(end)+taper]; + taper = taper*ratio; + end + if (stop - start_taper(end)) >= max_res + % fill with equidistant lines + n1 = ceil( (stop-start_taper(end)) / max_res ) + 1; + start_taper = [start_taper linspace(start_taper(end),stop,n1)]; + else + % not enough space for entire taper + start_taper(end) = []; % likely too near to stop + if numel(start_taper) > 0 + start_taper = [start_taper (stop+start_taper(end))/2]; % create a centered line + end + end + lines = sort(unique(start_taper(start_taper<=stop))); + return +end + +taper = start_res*ratio; +start_taper = start; +while (taper*ratio<max_res) && (start_taper(end)+taper < stop) + start_taper = [start_taper start_taper(end)+taper]; + taper = taper*ratio; +end +if (numel(start_taper) > 1) && (start_taper(end) - start_taper(end-1) > stop - start_taper(end)) + % not enough space for entire taper + start_taper(end) = []; % likely too near to stop + start_taper = [start_taper (stop+start_taper(end))/2]; % create a centered line +end + +taper = stop_res*ratio; +stop_taper = stop; +while (taper*ratio<max_res) && (stop_taper(1)-taper > start) + stop_taper = [stop_taper(1)-taper stop_taper]; + taper = taper*ratio; +end +if (numel(stop_taper) > 1) && (stop_taper(2) - stop_taper(1) > start - stop_taper(1)) + % not enough space for entire taper + stop_taper(1) = []; % likely too near to start + stop_taper = [(start+stop_taper(1))/2 stop_taper]; % create a centered line +end + +while ~isempty(start_taper) && ~isempty(stop_taper) && (stop_taper(1) < start_taper(end)) + + diff_start = diff(start_taper); + if isempty(diff_start) + diff_start = start_res; + end + diff_stop = diff(stop_taper); + if isempty(diff_stop) + diff_stop = stop_res; + end + + if (diff_start(end)>diff_stop(1)) + start_taper(end) = []; + else + stop_taper(1) = []; + end +end + +% it may happen, that the space between start_taper(end) and stop_taper(1) +% is very small; correct it: +if numel(start_taper)>=2 && numel(stop_taper)>=2 + d = diff( [start_taper(end-1) start_taper(end) stop_taper(1) stop_taper(2) ]); + if (d(1)/d(2) > ratio) || (d(3)/d(2) > ratio) + addLines = linspace( start_taper(end-1), stop_taper(2), 4 ); + lines = unique( [start_taper(1:end-1) addLines stop_taper(2:end)] ); + return + end +end + +if isempty(start_taper) + start_taper = start; +end +if isempty(stop_taper) + stop_taper = stop; +end + +% fill remaining space with equidistant lines +numL = ceil((stop_taper(1) - start_taper(end))/max_res)+1; +lines = unique([start_taper linspace(start_taper(end),stop_taper(1),numL) stop_taper]); +end % SmoothRange() |