function VisualizeSolution(U, PosX, PosY, MeshSize, t, Front, doVelocitiField, showBoundaryVelocities,QuiverSpec)   

UniqueX = unique(PosX);
UniqueY = unique(PosY);
NumNodesX = length(UniqueX);
NumNodesY = length(UniqueY);

[InsideMask, InteriorMask, BoundaryMask] = IdentifyNodeTypes(PosX, PosY, Front, NumNodesX);
InsideIdx = find(InsideMask);
BoundaryIdx = find(BoundaryMask);

if ~doVelocitiField
    U_local = U(InsideMask);
    if showBoundaryVelocities
        vel = zeros(size(Front.Polygon));
        for i = 1:size(Front.Polygon,1)
            vel(i,:) = FT_velo_func(t, Front.Polygon(i,:));
        end
    end
    
    tri = DelaunayTri([PosX(InsideMask)   PosY(InsideMask)]);
    [v, c] = voronoiDiagram(tri);
    
    cla;
    hold on
    tmp = zeros(1, NumNodesX*NumNodesY);
    tmp(InsideIdx) = 1:length(InsideIdx);
    I = tmp(InteriorMask);
    for i = I
        patch(v(c{i},1),v(c{i},2),U_local(i),'Edgecolor','none');
    end
    
    [NatNeighborBoundaryPairs, Vertices, VoronoiCells] =...
        GetVoronoiDecomposition(BoundaryIdx, PosX, PosY, InsideMask, [], NumNodesX, NumNodesY);
    [~, Directions, MidPoints] = GetBoundaryIndependentMetrics(NatNeighborBoundaryPairs, PosX, PosY);
    [~, ~, VerticesCoords,~] =...
        GetBoundaryVolumes(BoundaryIdx, MeshSize, VoronoiCells, Vertices, NatNeighborBoundaryPairs, Directions, MidPoints, Front,false);
    for j = 1:length(BoundaryIdx)
        Idx = BoundaryIdx(j);
        patch(VerticesCoords{j}(1,:),VerticesCoords{j}(2,:),U_local(InsideIdx == Idx),'Edgecolor','none');
    end
    plot([Front.Polygon(:,1); Front.Polygon(1,1)],[Front.Polygon(:,2); Front.Polygon(1,2)],'k-','LineWidth',1)
    
    if showBoundaryVelocities
        quiver(Front.Polygon(:,1),Front.Polygon(:,2),vel(:,1),vel(:,2))
        title(['V_m_a_x  = ' num2str(max(sqrt(sum(vel.^2,2))), 3)],'Fontsize', 14);
    end
else
    plot([Front.Polygon(:,1); Front.Polygon(1,1)],[Front.Polygon(:,2); Front.Polygon(1,2)],'k-','LineWidth',1)
    skip = round(NumNodesX/24);
    hold on
    i = 1:skip:NumNodesX;
    j = 1:skip:NumNodesY;
    I = repmat(i', 1, length(j));
    J = repmat(j, length(i),1);
    Idx = I(:) + NumNodesX*(J(:)-1);
    Mask = false(size(InsideMask));
    Mask(Idx) = true;
    Mask = Mask & InsideMask;
    Ux_local = U{1}(Mask);
    Uy_local = U{2}(Mask);
    quiver(PosX(Mask), PosY(Mask), Ux_local, Uy_local, U{3}, QuiverSpec{:})
    hold off
end
xlabel('X','Fontsize', 14);
ylabel('Y','Fontsize', 14);

hold off
