From dc411dfa4503c51caacaed4f2cfcdf2d06e41aa7 Mon Sep 17 00:00:00 2001
From: Ivar Stefansson <ist050@uib.no>
Date: Wed, 21 Dec 2016 18:29:20 +0100
Subject: [PATCH] Upload new file
---
scripts/compute_2d_error.m | 107 +++++++++++++++++++++++++++++++++++++
1 file changed, 107 insertions(+)
create mode 100644 scripts/compute_2d_error.m
diff --git a/scripts/compute_2d_error.m b/scripts/compute_2d_error.m
new file mode 100644
index 0000000..cf76305
--- /dev/null
+++ b/scripts/compute_2d_error.m
@@ -0,0 +1,107 @@
+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
+
--
GitLab