function[Data]=initializeInterface(Mesh, Data, it)

% initializeInterface: Initializes the interfaces properties
% INPUT:
% Mesh:     Structure containing all the mesh parameters
% Data:     Structure containing all the data parameters
% it:       Boundary tags for either a portion of the boundary or the entire boundary
%           The first component is the tag for body 1
%           The second component is the tag for body 2
% OUTPUT:
% Data:     Updated data structure with interface properties

% Retrieve mesh parameters
% Retrieve boundary tag
BC_tag=Mesh.BC_tag;
% Retrieve coordinate matrix
coord=Mesh.coord_mat;
% Define radius of intersection
if isfield(Data.Contact, 'radius')
    r=Data.Contact.radius;
else
    r=0.052;
end

% Boundary elements
elements1=Mesh.BC_nodes(ismember(BC_tag,it{1}),:);
elements2=Mesh.BC_nodes(ismember(BC_tag,it{2}),:);

% Boundary nodes
bc_nodes1=unique(elements1);
bc_nodes2=unique(elements2);

m1=length(bc_nodes1);
m2=length(bc_nodes2);

% Compute element lengths
l1=vecnorm(coord(elements1(:,2),:)-coord(elements1(:,1),:),2,2);
l2=vecnorm(coord(elements2(:,2),:)-coord(elements2(:,1),:),2,2);

% Compute mesh sizes
h1=min(l1);
h2=min(l2);

body_1.h=h1;
body_2.h=h2;

h=min([h1, h2]);
Data.Contact.h=h;

% Indicators
ic1=zeros(m1,1);
ic2=zeros(m2,1);

% Create interfaces based on closest nodes from opposing interface
for k=1:m1
    node=bc_nodes1(k);
    coord_node=coord(node,:);
    dist=vecnorm(coord(bc_nodes2,:)-coord_node,2,2);
    ic1(k)=any(dist<r);
    
end
for k=1:m2
    node=bc_nodes2(k);
    coord_node=coord(node,:);
    dist=vecnorm(coord(bc_nodes1,:)-coord_node,2,2);
    ic2(k)=any(dist<r);
end

ic1=logical(ic1);
ic2=logical(ic2);

interface1=bc_nodes1(ic1);
interface2=bc_nodes2(ic2);

% Correct interface nodes to account for finite elements
index1=all(ismember(Mesh.BC_nodes, interface1), 2);
index2=all(ismember(Mesh.BC_nodes, interface2), 2);
% Get interface connectivity matrices
body_1.interface_connect=Mesh.BC_nodes(index1,:);
body_2.interface_connect=Mesh.BC_nodes(index2,:);
% Get interface nodes
body_1.interface_nodes=unique(body_1.interface_connect);
body_2.interface_nodes=unique(body_2.interface_connect);
% Get interface degrees of freedom
interface_dof_1=Mesh.Dof_set(:,body_1.interface_nodes);
interface_dof_2=Mesh.Dof_set(:,body_2.interface_nodes);
body_1.interface_dof=interface_dof_1(:);
body_2.interface_dof=interface_dof_2(:);
% Initial active set
body_1.initial_nodes=body_1.interface_nodes;
body_2.initial_nodes=body_2.interface_nodes;
% Initial connectivity matrix
body_1.initial_interface_connect=body_1.interface_connect;
body_2.initial_interface_connect=body_2.interface_connect;
% Active set
body_1.active_set=ismember(body_1.initial_nodes, body_1.interface_nodes);
body_2.active_set=ismember(body_2.initial_nodes, body_2.interface_nodes);

% Initialization
body=cell(2,1);
% Create body entity
body{1}=body_1;
body{2}=body_2;
Data.body=body;
% Map global dof to local ones
[Data]=mapGlobal2Local(Data);
end