Newer
Older
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')
deltaP = max(p)-min(p);
if nargin == 2
kmat = 1;
end
fracIndices = find(k ~= kmat);
t(fracIndices) = [];
p(fracIndices) = [];
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
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