function [relativeMatrixError] = compute_2d_error( coarseFile, referenceFile, kmat )
%compute_2d_error Evaluates the norm in the entire matrix
%   Detailed explanation goes here
load(coarseFile,'T','P','Points')
Points = double(Points);
% Assumes the triangulation/connectivity list T is a matrix, i.e. all cells
% of equal number of vertices. Co-dimension one fracture cells may be left
% out (they shouldn't contribute to the norm anyway) if they cause trouble.
% Pressure vector P and npoints x 2 point list Points.

load(referenceFile,'t','p','x','k','domainarea')

points          = x;
deltaP = max(p)-min(p);
if nargin == 2
    kmat = 1;
end
fracIndices = find(k ~= kmat);
t(fracIndices) = [];
p(fracIndices) = [];
nCells          = length(t);
maxNNodes       = 4; % for the fine mesh. Reference solutions with both 
% quads and triangles make the code a bit messy.
tNan            = NaN(nCells,maxNNodes);
remainingCells  = true(nCells,1);
errorSquared    = zeros(nCells,1);

% Find the reference cell containing each of the coarse vertices
for c = 1:length(t)
    id                          = t{c};
    tNan(c,1:length(id))        = id;
    polygon                     = points(id,:);
    [in,on]     = inpolygon(Points(:,1),Points(:,2),polygon(:,1),polygon(:,2));
    PointId     = find(in & ~on);
    if ~isempty(PointId)
        if sum(PointId)>1
            UniquePoints = unique(1e-5.*round(1e5.*Points(PointId,:)),'rows'); %some methods 
            % define certain points twice
            if size(UniquePoints,1)>1
                disp(UniquePoints);
                errorSquared('small cell %i surrounds the above listed course points',c);
                % Occurs for too large fine cells compared to the coarse
                % cells. Violates assumptions used below and renders the
                % code useless.
            end
        end
        for PId = 1:length(PointId)
            Cells = find(any(T==PointId(PId),2));
            
            for i = 1:length(Cells)
                C       = Cells(i);
                Id      = T(C,:);
                Polygon = Points(Id,:);
                area    = area_on_point(polygon,Polygon,Points(PointId,:));
                dp      = P(C)-p(c);
                errorSquared(c)= errorSquared(c) + (dp)^2*area;
            end
        end
        remainingCells(c) = false;

    end
end
%%
% eliminate the cells already taken care of. 
t2              = t(remainingCells);
tNan2           = tNan(remainingCells,:);
nanInd          = isnan(tNan2(:,4)); % obs bare for 3+4
t3              = tNan2;
t3(nanInd,4)    = tNan2(nanInd,1);
errorSquared2   = zeros(size(tNan2,1),1);
p2              = p(remainingCells);


for C = 1:length(T)
    % construct coarse cell
    Id                  = T(C,:);
    Polygon             = Points(Id,:);

    % Find reference points inside the coarse cell
    [in,on]             = inpolygon(points(:,1),points(:,2),Polygon(:,1),Polygon(:,2));
    pointsIn            = find(in & ~on);

    % Distinguish between cells completely and partly inside large cell
    tIn                 = ismember(tNan2,pointsIn);
    cellsIn             = find(any(tIn,2));
    cellsCompletelyIn   = cellsIn(all(ismember(t3(cellsIn,:),pointsIn),2));
    cellsPartlyIn       = setdiff(cellsIn,cellsCompletelyIn);

    % Evaluate for the two cases
    for i = 1:length(cellsCompletelyIn)
        c               = cellsCompletelyIn(i);
        x               = points(t2{c},:);
        [~,area]        = convhull(x);
        dp              = (P(C)-p2(c))^2;
        errorSquared2(c)= dp*area;
    end
    for i = 1:length(cellsPartlyIn)
        c               = cellsPartlyIn(i);
        polygon         = points(t2{c},:);
        area            = area_on_face(polygon,Polygon);
        dp              = (P(C)-p2(c));
        errorSquared2(c)= errorSquared2(c) + dp^2*area;
    end
end

errorSquared(remainingCells) = errorSquared2;

relativeMatrixError = sqrt(sum(errorSquared))/(deltaP*sqrt(domainarea));



end