gridpart.m

function gridp = gridpart(grid,eind)
%function gridp = gridpart(grid,eind)
% function extracting a part of a .triagrid defined by the given 
% element indices in the vector 'eind'.
%
% @note The neighbour information of the new resulting boundaries is set to
% '-10'
%
% @note The properties gridbase.hmin, gridbase.alpha and the
% distance-information in the new boundary elements are simply copied. I.e.
% these fields do not completely meet the definition in the constructor. They
% might be chosen slightly different, such that the 'gridp' would be really
% identical to the result generated by the constructor on the subset of points.
%
% Parameters:
%  eind: vector of cell indices which shall be extracted from the grid.
%
% Return values:
%  gridp: the partial grid of type .gridbase with extracted cells 'eind'.

% Bernard Haasdonk 15.5.2007

gridp = triagrid(grid);
gridp.nelements = length(eind);

% generate elementid translation map
new_el_id = zeros(1,grid.nelements);
new_el_id(eind) = 1:length(eind);

% generate vertex translation map
new_vertex_id = zeros(1,grid.nvertices);
mask = zeros(1,grid.nvertices);
mask(grid.VI(eind,:))= 1;
vind = find(mask);
new_vertex_id(vind) = 1:length(vind);

gridp.nvertices = max(new_vertex_id);
gridp.A = gridp.A(eind);
gridp.Ainv = gridp.Ainv(eind);

gridp.X = grid.X(vind);
gridp.Y = grid.Y(vind);

gridp.VI = grid.VI(eind,:); %: VI(i,j) is the global vertex index of j-th
gridp.VI = new_vertex_id(gridp.VI);

gridp.CX = gridp.CX(eind);
gridp.CY = gridp.CY(eind);

gridp.NBI = grid.NBI(eind,:);
i = find(gridp.NBI>0);
gridp.NBI(i) = new_el_id(gridp.NBI(i));
i = find(gridp.NBI == 0);
if ~isempty(i)
  gridp.NBI(i)= -10;
end;

gridp.INB = grid.INB(eind,:);

gridp.EL = grid.EL(eind,:);
gridp.DC = grid.DC(eind,:);
gridp.NX = grid.NX(eind,:);
gridp.NY = grid.NY(eind,:);
gridp.ECX = grid.ECX(eind,:);
gridp.ECY = grid.ECY(eind,:);
gridp.SX = grid.SX(eind,:);
gridp.SY = grid.SY(eind,:);
gridp.ESX = grid.ESX(eind,:);
gridp.ESY = grid.ESY(eind,:);
gridp.DS = grid.DS(eind,:);
gridp.JIT = grid.JIT(eind,:,:);

gridp.nedges_boundary = length(find(gridp.NBI<0));
gridp.nedges_interior = 0.5*(3*gridp.nelements - gridp.nedges_boundary);