function [BoundaryIdxLocal] = GetBoundaryIndexLocal(InsideMask,size_U,BoundaryIdxGlobal)

  % returns the local boundary indices for the boundary nodes defined
  % globally in 'BoundaryIdxGlobal'. InsideMask is the mask corresponding
  % to inside nodes. 
  
  % First, creates a boundary mask 'BoundaryMask' with the same size as the
  % solution vector, filled with 'true' elements for boundary CVs.
  % 'size_U' has the size information of the full (includes inside and
  % outside nodes) variable.
  
  BoundaryMask = false(size_U);
  BoundaryMask(BoundaryIdxGlobal) = true;
  
  % Then create the corresponding index arrays for inside and boundary
  % nodes.
  InIndex = find(InsideMask);
  BIndex = find(BoundaryMask);

  % finally, create a local index of the boundary nodes
  BoundaryIdxLocal = ones(length(BIndex),1);
  for i=1:length(BIndex)
    BoundaryIdxLocal(i) = find(InIndex==BIndex(i));
  end
  
end