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 4495 2014-04-09 17:09:21Z 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 % Use a for loop, vectorised approach can run out of memory
0032 for i = 1:length(shapes)
0033    l1 = length(M1.nodes);
0034    M2 = shapes{i};
0035    if ~isfield(M2, 'boundary');
0036       M2.boundary = find_boundary(M2);
0037    end
0038    if ~isfield(M2, 'mat_idx')
0039       M2.mat_idx = {1:length(M2.elems)};
0040    end
0041    % make sure boundaries are the same dimension
0042    if size(M1.boundary,2) == 2 && size(M2.boundary,2)==3
0043       % M1 is 2D, M2 is 3D
0044       M1.boundary = M1.elems;
0045    elseif size(M1.boundary,2) == 3 && size(M2.boundary,2)==2
0046       M2.boundary = M2.elems;
0047    end
0048    
0049    useM1 = nodes_in_bounding_box(M2,M1,th);
0050    useM1 = find(useM1);
0051    useM2 = nodes_in_bounding_box(M1,M2,th);
0052    
0053    if isempty(useM1) || all(~useM2(:))
0054       M1 = naive_join(M1,M2);
0055       continue
0056    end
0057    
0058    nodes_to_add = M2.nodes(~useM2,:);
0059    n_new_nodes = nnz(~useM2);
0060    match = zeros(length(M2.nodes),1);
0061    match(~useM2) = l1 + (1:n_new_nodes);
0062    useM2 = find(useM2);
0063    for n = useM2'
0064       use = nodes_near_node(M1.nodes(useM1,:), M2.nodes(n,:), 1.1*th);
0065       switch nnz(use)
0066          case 0
0067          case 1
0068             match(n) = useM1(use);
0069          otherwise
0070             use = useM1(use);
0071             len_use = length(use);
0072             D = M1.nodes(use,:) - ones(len_use,1)*M2.nodes(n,:);
0073             D = sum(D.^2,2); % sqrt(sum(D.^2,2)); % don't need to sqrt
0074             [val p] = min(D);
0075             if val < th^2    % but then need th^2
0076                match(n) = use(p); % matching node of M1
0077             end
0078       end
0079       if ~match(n)
0080          match(n) = l1 + n_new_nodes+1;
0081          n_new_nodes = n_new_nodes + 1;
0082          nodes_to_add = [nodes_to_add; M2.nodes(n,:)];
0083       end
0084    end
0085    M1.nodes = [M1.nodes; nodes_to_add];
0086    LE = length(M1.elems);
0087    for j = 1:numel(M2.mat_idx)
0088       M1.mat_idx{end+1} = LE+M2.mat_idx{j};
0089    end
0090    %M1.mat_idx = [M1.mat_idx {length(M1.elems)+(1:length(M2.elems))}];
0091    M1.elems = [M1.elems; match(M2.elems)];
0092    % this is not strictly correct, but visualizes nicely
0093    M1.boundary = [M1.boundary; match(M2.boundary)]; 
0094 end
0095 
0096 out =  M1;
0097 % rmfield(M1,'boundary');
0098 
0099 function M1 = naive_join(M1,M2)
0100 LN = length(M1.nodes);
0101 LE = length(M1.elems);
0102 M1.nodes = [M1.nodes; M2.nodes];
0103 M1.elems = [M1.elems; LN+M2.elems];
0104 for i = 1:numel(M2.mat_idx)
0105    M1.mat_idx{end+1} = LE+M2.mat_idx{i};
0106 end
0107 M1.boundary = [M1.boundary; LN+M2.boundary];
0108    
0109 
0110 function use = nodes_in_bounding_box(M2,M1,th)
0111    % limit to nodes in M1 that are within the bounding box of M2
0112    maxM2 = max(M2.nodes)+th;
0113    minM2 = min(M2.nodes)-th;
0114    use = true(length(M1.nodes),1);
0115    for i = 1:size(M1.nodes,2);
0116       use = use & M1.nodes(:,i) < maxM2(i) & M1.nodes(:,i) > minM2(i);
0117    end
0118 
0119 function use = nodes_near_node(nodes,node,th)
0120    maxlim = node + th;
0121    minlim = node - th;
0122    use = true(length(nodes),1);
0123    for i = 1:size(nodes,2);
0124       use = use & nodes(:,i) < maxlim(i) & nodes(:,i) > minlim(i);
0125    end
0126 
0127 function run_unit_test
0128 subplot(221)
0129 cyl = ng_mk_cyl_models(3,[0],[]);
0130 show_fem(cyl)
0131 
0132 subplot(222)
0133 top_nodes = cyl.nodes(:,3)>=1.5;
0134 top_elems = sum(top_nodes(cyl.elems),2)==4;
0135 top.elems = cyl.elems(top_elems,:);
0136 nds = unique(top.elems);
0137 map = zeros(1,length(cyl.nodes));
0138 map(nds) = 1:length(nds);
0139 top.elems = map(top.elems);
0140 top.nodes = cyl.nodes(nds,:);
0141 top.type = 'fwd_model';
0142 top.boundary = find_boundary(top);
0143 show_fem(top)
0144 zlim([0 3]);
0145 
0146 subplot(223)
0147 bot_elems = ~top_elems;
0148 bot.elems = cyl.elems(bot_elems,:);
0149 nds = unique(bot.elems);
0150 map = zeros(1,length(cyl.nodes));
0151 map(nds) = 1:length(nds);
0152 bot.elems = map(bot.elems);
0153 bot.nodes = cyl.nodes(nds,:);
0154 bot.type = 'fwd_model';
0155 bot.boundary = find_boundary(bot);
0156 show_fem(bot)
0157 zlim([0 3]);
0158 
0159 
0160 subplot(224)
0161 M = merge_meshes(bot, top);
0162 show_fem(M);
0163 
0164 unit_test_cmp('Number of nodes',length(cyl.nodes), length(M.nodes),0);
0165 unit_test_cmp('Number of elems',length(cyl.elems), length(M.elems),0);