writeMesh/getInterf3d.m at main · mjx888/writeMesh · GitHub

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
function interfnode_cell = getInterf3d( nodecoor_cell )
% getInterf3d: get interface node between phase nodecoor_cell{i} and nodecoor_cell{j}
%
% interfnode_cell - n*n cell array. If nodecoor_cell{i} and
%                   nodecoor_cell{j} don't contact,  interfnode_cell{i,j}
%                   would be []. interfnode_cell{i,j} stores interfacial
%                   node number of nodecoor_cell{i} .
%
% Copyright (C) 2019-2025 by Jiexian Ma, mjx0799@gmail.com
% Distributed under the terms of the GNU General Public License (version 3)
%
% Project website: https://github.com/mjx888/im2mesh
%

    len = length( nodecoor_cell );
    interfnode_cell = cell( len, len );

    for i = 1: len-1
        for j = i+1: len
            nodecoor_A = nodecoor_cell{i};
            nodecoor_B = nodecoor_cell{j};
            [interfnode_A, interfnode_B] = findInterfNode(nodecoor_A, nodecoor_B);

            if isempty(interfnode_A) || isempty(interfnode_B)
                continue
            else
                interfnode_cell{i,j} = interfnode_A;
                interfnode_cell{j,i} = interfnode_B;
            end
        end
    end

end

function [interfnode_A, interfnode_B] = findInterfNode(nodecoor_A, nodecoor_B)
% find the interface between phase A and phase B, return node number in
% corresponding phase

    mark_A = false( length(nodecoor_A), 1 );
    mark_B = false( length(nodecoor_B), 1 );

    for i = 1: size( nodecoor_A( :, 2:4 ), 1 )
        idx = find( nodecoor_A( i, 2 ) == nodecoor_B( :, 2 ) & ...
                    nodecoor_A( i, 3 ) == nodecoor_B( :, 3 ) & ...
                    nodecoor_A( i, 4 ) == nodecoor_B( :, 4 ) );
        len = length(idx);

        switch len
            case 0
                continue
            case 1
                mark_A( i ) = 1;
                mark_B( idx ) =1;
            otherwise
                disp('multiple contact points')
        end
    end

    interfnode_A = nodecoor_A( mark_A, 1 ); % find 1
    interfnode_B = nodecoor_B( mark_B, 1 );

end