function u = solver(p,t,bpm,bnd,volumeforce,potential)

    coordinates = p; 
    elements3 = t; 
    elements4 = []; % quad nodes: we have none
    bpm.conditions  = []; %no special conditions for  bpm simulations

    % for simple shapes, set default boundary nodes as single vector
    if isempty(bnd.dirichlet) == 1, 
        bnd.dirichlet{1} = unique(boundedges(p,t));
    end

    % this vector contains ID of each dirichlet boundary.
    % we assume every boundary ad Dirichlet by default.
    dirichlet = 1:length(bnd.dirichlet); 
    neumann = 1:length(bnd.neumann); 

    % this vector contains id of each neumann boundary
    % and removes it from Dirichlet vector.
    %  for id = 1:length(neumann)
    %      dirichlet(find(dirichlet==neumann(id))) = []; 
    %  end

    % construct N and D closed boundaries by node pairs: 
    % [node1 node2; node2 node3; ... etc ... nodeN node1]
    dirichletpairs = [];
    neumannpairs = [];
    for nd = 1:length(dirichlet)
        dirichletpairs = [dirichletpairs; ConstructClosedBoundaryPairs(bnd.dirichlet{dirichlet(nd)})];
    end
    for nn = 1:length(neumann)
        neumannpairs = [neumannpairs; ConstructClosedBoundaryPairs(bnd.neumann{neumann(nn)})]; 
    end
    
    FreeNodes = setdiff(1:size(coordinates,1),unique(dirichletpairs));
    A = sparse(size(coordinates,1),size(coordinates,1));
    b = sparse(size(coordinates,1),1);

    % Assembly
    for j = 1:size(elements3,1)
        A(elements3(j,:),elements3(j,:)) = A(elements3(j,:), ...
        elements3(j,:)) + stima3(coordinates(elements3(j,:),:));
    end

    for j = 1:size(elements4,1)
        A(elements4(j,:),elements4(j,:)) = A(elements4(j,:), ...
        elements4(j,:)) + stima4(coordinates(elements4(j,:),:));
    end

    % Volume Forces (original)
    %for j = 1:size(elements3,1)
    %    b(elements3(j,:)) = b(elements3(j,:)) + ...
    %    det([1,1,1; coordinates(elements3(j,:),:)']) * ...
    %    f(VolumeForce, sum(coordinates(elements3(j,:),:))/3)/6;
    %end
    
    % Volume Forces (patched)
    for j = 1:size(elements3,1)
        centroid_coords   = sum(coordinates(elements3(j,:),:))/3;
        b(elements3(j,:)) = b(elements3(j,:)) + ...
                            det([1,1,1; coordinates(elements3(j,:),:)']) * ...
                            f(volumeforce, centroid_coords)/6 * rhs(centroid_coords,bpm.conditions);
    end

    for j = 1:size(elements4,1)
        b(elements4(j,:)) = b(elements4(j,:)) + ...
        det([1,1,1; coordinates(elements4(j,1:3),:)']) * ...
        f(volumeforce, sum(coordinates(elements4(j,:),:))/4)/4;
    end

    % Neumann conditions
    for j = 1 : size(neumannpairs,1)
        b(neumannpairs(j,:))=b(neumannpairs(j,:)) + ...
        norm(coordinates(neumannpairs(j,1),:)-coordinates(neumannpairs(j,2),:))*...
        g(sum(coordinates(neumannpairs(j,:),:))/2)/2;
    end

    % Dirichlet conditions
    u = sparse(size(coordinates,1),1);
    % u(unique(dirichletpairs)) = u_d(coordinates(unique(dirichletpairs),:)); 
    for nd = 1:length(dirichlet)
        u(bnd.dirichlet{dirichlet(nd)}) = u_d_v(p(bnd.dirichlet{dirichlet(nd)},:),potential(dirichlet(nd)));
    end 
    b = b - A * u;

    % Computation of the solution
    u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);

    % graphic representation
    %show(elements3,elements4,coordinates,full(u),geom);
    
% Some other auxiliar functions

function M = stima3(vertices)
    d = size(vertices,2);
    G = [ones(1,d+1); vertices'] \ [zeros(1,d); eye(d)];
    M = det([ones(1,d+1); vertices']) * (G * G') / prod(1:d);
end

function M = stima4(vertices)
    D_Phi = [vertices(2,:)-vertices(1,:); vertices(4,:) - vertices(1,:)]';
    B = inv(D_Phi' * D_Phi);
    C1 = [2,-2;-2,2]*B(1,1) + [3,0;0,-3]*B(1,2) + [2,1;1,2]*B(2,2);
    C2 = [-1,1;1,-1]*B(1,1) + [-3,0;0,3]*B(1,2) + [-1,-2;-2,-1]*B(2,2);
    M = det(D_Phi) * [C1 C2; C2 C1] / 6;
end

function DirichletBoundaryValue = u_d(x)
    DirichletBoundaryValue = zeros(size(x,1),1);
end

function DirichletBoundaryValue = u_d_v(x,v)
    % modified u_d which accepts potential value for each boundary
    DirichletBoundaryValue = zeros(size(x,1),1) + v;
end

function VolumeForce = f(VolumeForce,x)
    VolumeForce = VolumeForce + zeros(size(x,1),1);
end

function bla = rhs(centroid_coords, conditions)
    if strcmp(conditions,'square_space_charge') == 1
        bla = 2*centroid_coords(1)*(1-centroid_coords(1))...
             +2*centroid_coords(2)*(1-centroid_coords(2));
    else bla = 1;
    end
            
end

function Stress = g(x)
    Stress = zeros(size(x,1),1);
end
    
function b = ConstructClosedBoundaryPairs(a)
    for i = 1:length(a), 
        if i == length(a), 
            b(i,:) = [a(end) a(1)]; % close the polygon: last + first
        else
            b(i,:) = [a(i) a(i+1)]; % build polygon edges
        end 
    end 
end

end
