function [A,BoundaryIdxLocal] = ApplyDirichletBndCondMat(A,InsideMask,NumInsideNodes,size_U,BoundaryIdxGlobal)

  % applies Dirichlet boundary condition to a matrix A for the boundary
  % nodes defined in 'BoundaryIdxGlobal'. InsideMask is the mask
  % corresponding to inside nodes. NumInsideNodes is the size of the
  % corresponding system (inside nodes only). size_U is the size of the
  % corresponding state variable vector (includes both inside and outside
  % nodes). 
  %
  % First, create a boundary mask 'BoundaryMask' with the same size as the
  % solution vector 'size_U', and fill it with 'true' elements for boundary
  % CVs.  
  BoundaryMask = false(size_U);
  BoundaryMask(BoundaryIdxGlobal) = true;
  
  % Then create the corresponding index arrays for inside and boundary
  % nodes.
  InIndex = find(InsideMask);
  BIndex = find(BoundaryMask);

  % 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
  
  % final matrix manipulation
  for i=1:length(BoundaryIdxLocal)
      A(BoundaryIdxLocal(i),1:NumInsideNodes) = 0;
      A(BoundaryIdxLocal(i),BoundaryIdxLocal(i)) = 1;
  end    

end