merge_meshes

MERGE_MESHES - merges two meshes sharing only boundary 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 (sequentially)

 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 sharing only boundary nodes
0003 % MERGE_MESHES(M1,M2,T) merges M2 in M1 using threshold T for
0004 %     detecting corresponding nodes. The meshes must not overlap.
0005 % MERGE_MESHES(M1,M2,M3,..., T) merges M2, M3, ... into M1 (sequentially)
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, 2013-2024. Licence: GPL v2 or v3
0014 % $Id: merge_meshes.m 6984 2024-11-14 20:24:35Z 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 
0026 if ~isfield(M1, 'mat_idx') || isempty(M1.mat_idx)
0027    idx = 1:uint32(length(M1.elems));
0028    M1.mat_idx = {idx(:)};
0029 end
0030 
0031 
0032 out = M1;
0033 out.elems = uint32(out.elems);
0034 try out.boundary = uint32(out.boundary); end
0035 if ~isfield(out, 'boundary') || isempty(out.boundary)
0036    out.boundary = find_boundary(out);
0037 end
0038    
0039 
0040 for i = 1:length(shapes)
0041     if length(shapes) > 1
0042        msg = sprintf('Merging mesh %d/%d ... ',i,length(shapes));
0043     else
0044        msg = 'Merging meshes ... ';
0045     end
0046     progress_msg(msg);
0047     
0048     M2 = shapes{i};
0049     M2.elems = uint32(M2.elems);
0050     try M2.boundary = uint32(M2.boundary); end
0051 
0052     if ~isfield(M2, 'boundary') || isempty(M2.boundary)
0053         M2.boundary = find_boundary(M2);
0054     end
0055     if ~isfield(M2, 'mat_idx') || isempty(M2.mat_idx)
0056         M2.mat_idx = {1:uint32(length(M2.elems))};
0057     end
0058     
0059     if size(out.elems,2) == 3 % surface mesh
0060         B1 = unique(out.elems);
0061     else
0062         B1 = unique(find_boundary(out));
0063     end
0064     
0065     if size(M2.elems,2) == 3 %surface mesh
0066         B2 = unique(M2.elems);
0067     else
0068         B2 = unique(find_boundary(M2));
0069     end
0070     
0071     M1_num_nodes = uint32(size(out.nodes,1));
0072     M2_num_nodes = uint32(size(M2.nodes,1));
0073     
0074     excl_B1 = ~nodes_in_bounding_box(M2.nodes(B2,:), out.nodes(B1,:), th);
0075     excl_B2 = ~nodes_in_bounding_box(out.nodes(B1,:), M2.nodes(B2,:), th);
0076     
0077     if nnz(excl_B1) == numel(B1) || nnz(excl_B2) == numel(B2)
0078         out.nodes = [out.nodes; M2.nodes];
0079         nodemap = M1_num_nodes + (1:M2_num_nodes);
0080     else
0081         B1(excl_B1) = [];
0082         B2(excl_B2) = [];
0083         
0084         bnd_nodes = false(size(M2.nodes,1),1);
0085         bnd_nodes(B2) = true;
0086         
0087         % add not-B2 to M1
0088         out.nodes = [out.nodes; M2.nodes(~bnd_nodes,:)];
0089         % build nodemap for B2
0090         nodemap = zeros(M2_num_nodes,1,'uint32');
0091         nodemap(~bnd_nodes) = M1_num_nodes + (1:uint32(M2_num_nodes - length(B2)));
0092         % find closest points in B1
0093         
0094         mem_use = 256 * 1024^2; %MiB
0095         if ispc || exist('OCTAVE_VERSION','builtin') &&  isunix
0096            meminfo = memory;
0097            mem_use = meminfo.MaxPossibleArrayBytes / 4;
0098         end
0099         % compute distance matrix in chunks to prevent exhausting memory
0100         chunk_sz = floor(mem_use / 8 / numel(B1)); % mem_use / 8 bytes(double)
0101         n_chunks = ceil(numel(B2)/chunk_sz);
0102         m = zeros(1,numel(B2));
0103         p = zeros(1,numel(B2));
0104         
0105         for c = 1:n_chunks
0106             progress_msg(c/n_chunks);
0107             idx = ((c-1)*chunk_sz + 1) : min(c*chunk_sz, numel(B2));
0108             D = (out.nodes(B1,1) - M2.nodes(B2(idx),1)').^2;
0109             for dim = 2:size(out.nodes,2)
0110                 D = D + (out.nodes(B1,dim) - M2.nodes(B2(idx),dim)').^2;
0111             end
0112             [m(idx), p(idx)] = min(D,[],1);
0113         end
0114         
0115         new_bnd_nodes = m > th^2;
0116         out.nodes = [out.nodes; M2.nodes(B2(new_bnd_nodes),:)];
0117         nodemap(B2(new_bnd_nodes)) = max(M1_num_nodes, max(nodemap)) + (1:uint32(nnz(new_bnd_nodes)));
0118         
0119         B2(new_bnd_nodes) = [];
0120         p(new_bnd_nodes) = [];
0121         
0122         nodemap(B2) = B1(p); 
0123     end
0124     
0125     M2.elems    = nodemap(M2.elems);
0126     M2.boundary = nodemap(M2.boundary);
0127     if isfield(M2, 'electrode')
0128         for e = 1:numel(M2.electrode)
0129             try
0130             M2.electrode(e).nodes = nodemap(M2.electrode(e).nodes);
0131             end
0132             try
0133             M2.electrode(e).faces = nodemap(M2.electrode(e).faces);
0134             end
0135         end
0136     end
0137     LE = length(out.elems);
0138     out.elems = [out.elems; M2.elems];
0139     d1 = size(out.boundary, 2);
0140     d2 = size(M2.boundary, 2);
0141     if d1 < d2 % out is 2D or surface mesh
0142        out.boundary = out.elems;
0143     end
0144     if d2 < d1
0145        M2.boundary = M2.elems;
0146     end
0147     out.boundary = unique([out.boundary; M2.boundary], 'rows');
0148     if isfield(M2, 'electrode') && isstruct(M2.electrode)
0149         if ~isfield(out, 'electrode'), out.electrode = []; end
0150         n = numel(out.electrode);
0151         % is there a better way to do this?
0152         fn = fieldnames(M2.electrode);
0153         for j = 1:numel(M2.electrode)
0154             for f = 1:numel(fn)
0155                 out.electrode(n+j).(fn{f}) =  M2.electrode(j).(fn{f});
0156             end
0157         end
0158     end
0159     for j = 1:numel(M2.mat_idx)
0160         out.mat_idx{end+1} = LE+M2.mat_idx{j};
0161     end
0162     progress_msg(Inf);
0163 end
0164 
0165 % standard field order
0166 out = eidors_obj('fwd_model', out);
0167 
0168 function use = nodes_in_bounding_box(LIMnodes,Mnodes,th)
0169    % limit to nodes in M1 that are within the bounding box of M2
0170    maxLIM = max(LIMnodes)+th;
0171    minLIM = min(LIMnodes)-th;
0172    use = true(length(Mnodes),1);
0173    for i = 1:size(Mnodes,2)
0174       use = use & Mnodes(:,i) < maxLIM(i) & Mnodes(:,i) > minLIM(i);
0175    end
0176 
0177 function run_unit_test
0178 subplot(221)
0179 cyl = ng_mk_cyl_models(3,[0],[]);
0180 show_fem(cyl)
0181 
0182 subplot(222)
0183 top_nodes = cyl.nodes(:,3)>=1.5;
0184 top_elems = sum(top_nodes(cyl.elems),2)==4;
0185 top.elems = cyl.elems(top_elems,:);
0186 nds = unique(top.elems);
0187 map = zeros(1,length(cyl.nodes));
0188 map(nds) = 1:length(nds);
0189 top.elems = map(top.elems);
0190 top.nodes = cyl.nodes(nds,:);
0191 top.type = 'fwd_model';
0192 top.boundary = find_boundary(top);
0193 show_fem(top)
0194 zlim([0 3]);
0195 
0196 subplot(223)
0197 bot_elems = ~top_elems;
0198 bot.elems = cyl.elems(bot_elems,:);
0199 nds = unique(bot.elems);
0200 map = zeros(1,length(cyl.nodes));
0201 map(nds) = 1:length(nds);
0202 bot.elems = map(bot.elems);
0203 bot.nodes = cyl.nodes(nds,:);
0204 bot.type = 'fwd_model';
0205 bot.boundary = find_boundary(bot);
0206 show_fem(bot)
0207 zlim([0 3]);
0208 
0209 
0210 subplot(224)
0211 M = merge_meshes(bot, top);
0212 show_fem(M);
0213 
0214 unit_test_cmp('Number of nodes',length(cyl.nodes), length(M.nodes),0);
0215 unit_test_cmp('Number of elems',length(cyl.elems), length(M.elems),0);