function [u,K,M,b]=FEPoissonPenalized(f,g,N,T);
% FEPOISSON solves the Poisson equation using finite elements
% with penalized Dirichlet condition
% u=FEPoisson(f,g,N,T); solves the Poisson equation on the
% triangular mesh given by the list of triangles T and nodes N
% with Dirichlet boundary conditions given by the function g and
% forcing function f.
n=size(T,1); m=size(N,1);
bn=zeros(m,1);
K=sparse(m,m); M=sparse(m,m);
tgv=10^(20);
%tgv=10^(10);

for i=1:n,
    Ke=ComputeElementStiffnessMatrix([N(T(i,1),:); N(T(i,2),:); N(T(i,3),:)]);
    Me=ComputeElementMassMatrix([N(T(i,1),:); N(T(i,2),:); N(T(i,3),:)]);
    bn(T(i,1))=bn(T(i,1)) | T(i,4) | T(i,6); % on the boundary
    bn(T(i,2))=bn(T(i,2)) | T(i,4) | T(i,5);
    bn(T(i,3))=bn(T(i,3)) | T(i,5) | T(i,6);
    K(T(i,1:3),T(i,1:3))=K(T(i,1:3),T(i,1:3))+Ke; % assemble
    M(T(i,1:3),T(i,1:3))=M(T(i,1:3),T(i,1:3))+Me;
end;
b=feval(f,N(:,1),N(:,2)); % right hand side and
b=M*b;
% if ind is a table of type logical (with   0 and  1), 
%then K(ind, :) keeps only the lines j with ind(j) = 1.
% and it does it in sparse
bl=logical(bn);
bln=~logical(bn);
Kd=K(bl,bl);
K(bl,bl)=Kd+ tgv*eye(size(Kd));

%gd=zeros(m,1);
% evaluate the boundary condition only on the boundary nodes
gd=feval(g,N(bl,1),N(bl,2));
b(bl)=b(bl)+tgv*gd;
warning('off','MATLAB:nearlySingularMatrix')

u=K\b;
warning('on','MATLAB:nearlySingularMatrix')
