Upload new file

scripts/compute_2d_error.m

+function [relativeMatrixError] = compute_2d_error( coarseFile,referenceFile )
+%compute1dError 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','domainarea')
+
+points          = x;
+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;
+
+deltaP = max(p)-min(p);
+relativeMatrixError = sqrt(sum(errorSquared))/(deltaP*sqrt(domainarea));
+
+
+
+end
+