merge_meshes

MERGE_MESHES - merges two meshes based on common nodes
 MERGE_MESHES(M1,M2,T) merges M2 in M1 using threshold T for detecting
     corresponding nodes. The meshes must not overlap.
 MERGE_MESHES(M1,M2,M3,..., T) merges M2, M3, ... into M1 (individually)

 Note that the boundaries of the separate meshes will only be
 concatenated, as this visualises nicely. To calculate the correct
 boundary use FIND_BOUNDARY.

 See also FIND_BOUNDARY

0001 function out = merge_meshes(M1,varargin)
0002 %MERGE_MESHES - merges two meshes based on common nodes
0003 % MERGE_MESHES(M1,M2,T) merges M2 in M1 using threshold T for detecting
0004 %     corresponding nodes. The meshes must not overlap.
0005 % MERGE_MESHES(M1,M2,M3,..., T) merges M2, M3, ... into M1 (individually)
0006 %
0007 % Note that the boundaries of the separate meshes will only be
0008 % concatenated, as this visualises nicely. To calculate the correct
0009 % boundary use FIND_BOUNDARY.
0010 %
0011 % See also FIND_BOUNDARY
0012 
0013 % (C) Bartlomiej Grychtol and Andy Adler, 2012-2013. Licence: GPL v2 or v3
0014 % $Id: merge_meshes.m 6454 2022-12-04 22:11:44Z bgrychtol $
0015 
0016 if ischar(M1) && strcmp(M1,'UNIT_TEST'), run_unit_test; return; end
0017 
0018 if nargin < 3  || isstruct(varargin{end})
0019    th = mean(std(M1.nodes))/length(M1.nodes);
0020    shapes = varargin;
0021 else
0022    th = varargin{end};
0023    shapes = varargin(1:end-1);
0024 end
0025 if ~isfield(M1, 'boundary');
0026    M1.boundary = find_boundary(M1);
0027 end
0028 if ~isfield(M1, 'mat_idx')
0029    M1.mat_idx = {1:length(M1.elems)};
0030 end
0031 
0032 do_msg = length(shapes) > 1;
0033 if do_msg, progress_msg(sprintf('Merging mesh %d/%d ... ',0,length(shapes))); end
0034 
0035 % Use a for loop, vectorised approach can run out of memory
0036 for i = 1:length(shapes)
0037    if do_msg
0038       progress_msg(sprintf('Merging mesh %d/%d ... ',i,length(shapes)));
0039    end
0040    l1 = length(M1.nodes);
0041    M2 = shapes{i};
0042    if ~isfield(M2, 'boundary');
0043       M2.boundary = find_boundary(M2);
0044    end
0045    if ~isfield(M2, 'mat_idx')
0046       M2.mat_idx = {1:length(M2.elems)};
0047    end
0048    % make sure boundaries are the same dimension
0049    if size(M1.boundary,2) == 2 && size(M2.boundary,2)==3
0050       % M1 is 2D, M2 is 3D
0051       M1.boundary = M1.elems;
0052    elseif size(M1.boundary,2) == 3 && size(M2.boundary,2)==2
0053       M2.boundary = M2.elems;
0054    end
0055    
0056    useM1 = nodes_in_bounding_box(M2,M1,th);
0057    useM1 = find(useM1);
0058    useM2 = nodes_in_bounding_box(M1,M2,th);
0059    
0060    if isempty(useM1) || all(~useM2(:))
0061       M1 = naive_join(M1,M2);
0062       continue
0063    end
0064    
0065    nodes_to_add = M2.nodes(~useM2,:);
0066    n_new_nodes = nnz(~useM2);
0067    match = zeros(length(M2.nodes),1);
0068    match(~useM2) = l1 + (1:n_new_nodes);
0069    useM2 = find(useM2);
0070    N = length(useM2');
0071    for j = 1:N
0072       if mod(j,100)==1, progress_msg(j/N); end
0073       n = useM2(j); 
0074       
0075       use = nodes_near_node(M1.nodes(useM1,:), M2.nodes(n,:), 1.1*th);
0076       switch nnz(use)
0077          case 0
0078          case 1
0079             match(n) = useM1(use);
0080          otherwise
0081             use = useM1(use);
0082             len_use = length(use);
0083             D = M1.nodes(use,:) - ones(len_use,1)*M2.nodes(n,:);
0084             D = sum(D.^2,2); % sqrt(sum(D.^2,2)); % don't need to sqrt
0085             [val p] = min(D);
0086             if val < th^2    % but then need th^2
0087                match(n) = use(p); % matching node of M1
0088             end
0089       end
0090       if ~match(n)
0091          match(n) = l1 + n_new_nodes+1;
0092          n_new_nodes = n_new_nodes + 1;
0093          nodes_to_add = [nodes_to_add; M2.nodes(n,:)];
0094       end
0095    end
0096    if do_msg, progress_msg(Inf); end
0097    M1.nodes = [M1.nodes; nodes_to_add];
0098    LE = length(M1.elems);
0099    for j = 1:numel(M2.mat_idx)
0100       M1.mat_idx{end+1} = LE+M2.mat_idx{j};
0101    end
0102    %M1.mat_idx = [M1.mat_idx {length(M1.elems)+(1:length(M2.elems))}];
0103    M1.elems = [M1.elems; match(M2.elems)];
0104    % this is not strictly correct, but visualizes nicely
0105    M1.boundary = [M1.boundary; match(M2.boundary)]; 
0106 end
0107 
0108 out =  M1;
0109 % rmfield(M1,'boundary');
0110 
0111 function M1 = naive_join(M1,M2)
0112 LN = length(M1.nodes);
0113 LE = length(M1.elems);
0114 M1.nodes = [M1.nodes; M2.nodes];
0115 M1.elems = [M1.elems; LN+M2.elems];
0116 for i = 1:numel(M2.mat_idx)
0117    M1.mat_idx{end+1} = LE+M2.mat_idx{i};
0118 end
0119 M1.boundary = [M1.boundary; LN+M2.boundary];
0120    
0121 
0122 function use = nodes_in_bounding_box(M2,M1,th)
0123    % limit to nodes in M1 that are within the bounding box of M2
0124    maxM2 = max(M2.nodes)+th;
0125    minM2 = min(M2.nodes)-th;
0126    use = true(length(M1.nodes),1);
0127    for i = 1:size(M1.nodes,2);
0128       use = use & M1.nodes(:,i) < maxM2(i) & M1.nodes(:,i) > minM2(i);
0129    end
0130 
0131 function use = nodes_near_node(nodes,node,th)
0132    maxlim = node + th;
0133    minlim = node - th;
0134    use = true(size(nodes,1),1);
0135    for i = 1:size(nodes,2);
0136       use = use & nodes(:,i) < maxlim(i) & nodes(:,i) > minlim(i);
0137    end
0138 
0139 function run_unit_test
0140 subplot(221)
0141 cyl = ng_mk_cyl_models(3,[0],[]);
0142 show_fem(cyl)
0143 
0144 subplot(222)
0145 top_nodes = cyl.nodes(:,3)>=1.5;
0146 top_elems = sum(top_nodes(cyl.elems),2)==4;
0147 top.elems = cyl.elems(top_elems,:);
0148 nds = unique(top.elems);
0149 map = zeros(1,length(cyl.nodes));
0150 map(nds) = 1:length(nds);
0151 top.elems = map(top.elems);
0152 top.nodes = cyl.nodes(nds,:);
0153 top.type = 'fwd_model';
0154 top.boundary = find_boundary(top);
0155 show_fem(top)
0156 zlim([0 3]);
0157 
0158 subplot(223)
0159 bot_elems = ~top_elems;
0160 bot.elems = cyl.elems(bot_elems,:);
0161 nds = unique(bot.elems);
0162 map = zeros(1,length(cyl.nodes));
0163 map(nds) = 1:length(nds);
0164 bot.elems = map(bot.elems);
0165 bot.nodes = cyl.nodes(nds,:);
0166 bot.type = 'fwd_model';
0167 bot.boundary = find_boundary(bot);
0168 show_fem(bot)
0169 zlim([0 3]);
0170 
0171 
0172 subplot(224)
0173 M = merge_meshes(bot, top);
0174 show_fem(M);
0175 
0176 unit_test_cmp('Number of nodes',length(cyl.nodes), length(M.nodes),0);
0177 unit_test_cmp('Number of elems',length(cyl.elems), length(M.elems),0);