% separate the boundary in Dirichlet and Neumann parts

D_index = find(IX.coord(:,1)<=x_level); % Dirichlet
N_index = find(IX.coord(:,1)> x_level); % Neumann

% Dirichlet part                  % Neumann b.c. part
IX_D.coord = IX.coord(D_index,:); IX_N.coord = IX.coord(N_index,:);
IX_D.type  = IX.type(D_index,:);  IX_N.type  = IX.type(N_index,:); 
IX_D.anch  = IX.anch(D_index,:);  IX_N.anch  = IX.anch(N_index,:); 
IX_D.n     = IX.n(D_index,:);     IX_N.n     = IX.n(N_index,:); 
IX_D.t     = IX.t(D_index,:);     IX_N.t     = IX.t(N_index,:);

% reassign the order:
IX.coord = [IX_D.coord; IX_N.coord];
IX.type  = [IX_D.type ; IX_N.type ];
IX.anch  = [IX_D.anch ; IX_N.anch ];
IX.n     = [IX_D.n    ; IX_N.n    ];
IX.t     = [IX_D.t    ; IX_N.t    ];
