summaryrefslogtreecommitdiff
path: root/CSXCAD/matlab/RecursiveSmoothMesh.m
blob: 0ce791c0f28142d67e5341d68e5ddf4aa37b7a58 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
function resultMeshPoints = RecursiveSmoothMesh(fixedMeshPoints, max_res, ratio, varargin)
% function resultMeshPoints = RecursiveSmoothMesh(fixedMeshPoints, max_res, ratio, varargin)
%
% RecursiveSmoothMesh: This function tries to generate an optimal mesh
%  between the given fixedMeshPoints. The space between the fixed points is
%  filled with the largest possible mesh step size considdering mesh max
%  resolution and mesh increase/decrease ratio.
%
% Algorithm: Mesh creation is done via a simple try and error approach:
%  The next mesh point could be between lastDistance*ratio and
%  lastDistance/ratio away from the last mesh point. If a fixed mesh point
%  exist in this distance, this point is used. If there is a fixed mesh
%  point below this distance, the last (determined) mesh point is wrong, go
%  back and try with a mesh point a few mm below. (If that is not possible
%  because the distance to the mesh point behind the last point is to low,
%  also update this mesh point).
%  The algorithm seems to work very well, except if the ratio is to low.
%  At a value of 1.2, the calculation takes a very long time.
%
%
% Parameter:    fixedMeshPoints:    List containing points at which a mesh
%                                    line should appear
%               max_res:            Maximum distance between two mesh lines
%               ratio:              Maximum grading ratio between two
%                                    neighbour mesh lines
%
% 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)
%
% Returns:      resultMeshPoints:   List containing the positions of all
%                                     mesh lines. If a empty list is
%                                     returned, no resolution could be
%                                     found.
%
% CSXCAD matlab interface
% -----------------------
% Author: Florian Pfanner
% Date: 28.09.2012


if (nargin < 3)
    ratio = 1.4;
end

if (ratio < 1.2)
    warning('ratio must be > 1.2, set it to 1.2');
    ratio = 1.2;
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

fixedMeshPoints = sort(unique(fixedMeshPoints(:)'));
if (length(fixedMeshPoints) < 2)
    error('Not enoughts points to create mesh!');
end


% special behaviour to generate a symetric mesh (If the fixedMeshPoints are
% given in a symetric way.)
symMesh = CheckSymmtricLines(fixedMeshPoints);
if (symMesh)
    nPoints = length(fixedMeshPoints);

    symMeshEven = mod(nPoints,2)==0;

    if (symMeshEven)
        center = mean(fixedMeshPoints);
        if (center-fixedMeshPoints(nPoints/2) > 0.5*max_res)
            fixedMeshPoints = [fixedMeshPoints(1:nPoints/2) center];
        else
            % this is not working!!!
            % fixedMeshPoints = fixedMeshPoints(1:nPoints/2+1);
            fixedMeshPoints = [fixedMeshPoints(1:nPoints/2) center];
        end
    else
        % For symetric mesh, only the first half of the fixed mesh points is
        % needed. But add center of the symetric mesh to the fixed mesh points
        % to ensure the center is included.
        fixedMeshPoints = fixedMeshPoints(1:(nPoints+1)/2);
    end
end


minDistance = min(fixedMeshPoints(2:end)-fixedMeshPoints(1:end-1));

% two cases:
%  * minDistance is smaller than max_res
%       -> try spaces for the first mesh between max_res and minDistance/2
%  * max_res is smaller than minDistance
%       -> try spaces for the first mesh between max_res and max_res/10
if (minDistance < max_res)
    trySpaces = linspace(max_res, minDistance/2, 10);
else
    trySpaces = linspace(max_res, max_res/10, 10);
end
for k=1:length(trySpaces)
    [resultMeshPoints, success] = recursiveMeshSearch(fixedMeshPoints(1), ...
        trySpaces(k), fixedMeshPoints(2:end), max_res, ratio);
    if (success)
        resultMeshPoints = [fixedMeshPoints(1) resultMeshPoints];
        break;
    end
end


if (symMesh)
    resultMeshPoints = [resultMeshPoints(1:end-1) max(resultMeshPoints)*2-resultMeshPoints(end:-1:1)];
end

% check result
if (check_mesh)
    CheckMesh(resultMeshPoints,0,max_res,ratio,0);
end

end % main function

function [meshPoints, success] = recursiveMeshSearch(lastMeshPoint, ...
    lastSpace, fixedMeshPoints, max_res, ratio)

nextMeshPointMax = lastMeshPoint + min(lastSpace*ratio, max_res);
nextMeshPointMin = lastMeshPoint + lastSpace/ratio;

% check if below ratio is a fixed mesh point -> abbort
if (fixedMeshPoints(1) < nextMeshPointMin)
    meshPoints = [];
    success = false;
    return;
end

% check if in range of ratio is a fixed mesh point:
if (fixedMeshPoints(1) >= nextMeshPointMin && fixedMeshPoints(1) <= nextMeshPointMax)
    % yes, use this fixed mesh point
    nextMeshPoint = fixedMeshPoints(1);
    if (length(fixedMeshPoints) > 1)
        [meshPointsCall, success] = recursiveMeshSearch(nextMeshPoint, ...
            nextMeshPoint-lastMeshPoint, fixedMeshPoints(2:end), ...
            max_res, ratio);
        if (success)    % backtracking was successful, return result
            meshPoints = [nextMeshPoint meshPointsCall];
        else            % no possible resulution found, last step has to be repeated!
            meshPoints = [];
        end
    else            % this was the last mesh point, finish!
        meshPoints = nextMeshPoint;
        success = true;
    end
    return;
end

% try to set next mesh point, begin with highest width
trySpace = linspace(nextMeshPointMax, nextMeshPointMin, 10);
for i=1:length(trySpace)
    [meshPointsCall, success] = recursiveMeshSearch(trySpace(i), ...
        trySpace(i)-lastMeshPoint, fixedMeshPoints, max_res, ...
        ratio);
    if (success)    % backtracking was successful, return
        meshPoints = [trySpace(i) meshPointsCall];
        return;
    end
end
% no sucessful points found, return false
meshPoints = [];
success = false;
end