[scripts] calculate the error only on the exact matrix domain

scripts/compute_1d_error.m

@@ -7,8 +7,13 @@ function [totalRelativeFractureError] = compute_1d_error( coarseFile,referenceFi
 % The reference file should contain a triangulation t, point list x and
 % solution vector p.
 
-load(referenceFile,'p','t','x')
+load(referenceFile,'p','t','x','k')
 load(coarseFile,'fracXAndPCell')
+maxp = max(p);
+minp = min(p);
+matrixIndices = find(k == 1);
+t(matrixIndices) = [];
+p(matrixIndices) = [];
 
 nFractures                      = length(fracXAndPCell);
 idividualRelativeFractureError  = zeros(nFractures,1); %Included for comparison
@@ -20,7 +25,7 @@ frac = true;
 for i = 1:nFractures
     fracXAndP                   = fracXAndPCell{i};
     [idividualRelativeFractureError(i),fractureErrorSquared(i),normalizationSquared(i)] = ...
-                                  line_norm(fracXAndP,p,t,x,frac);
+                                  line_norm(fracXAndP,p,t,x,frac,maxp,minp);
 end
 
 totalRelativeFractureError = sqrt(sum(fractureErrorSquared)/sum(normalizationSquared));

scripts/compute_2d_error.m

@@ -8,9 +8,13 @@ Points = double(Points);
 % 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')
+load(referenceFile,'t','p','x','k','domainarea')
 
 points          = x;
+deltaP = max(p)-min(p);
+fracIndices = find(k ~= 1);
+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.
@@ -98,7 +102,6 @@ end
 
 errorSquared(remainingCells) = errorSquared2;
 
-deltaP = max(p)-min(p);
 relativeMatrixError = sqrt(sum(errorSquared))/(deltaP*sqrt(domainarea));
 
 

scripts/hdf5_to_mat_mixed.m

+function hdf5_to_mat_matrix(infilename, cellwise, solutionname, outfilename)
+% hdf5_to_mat_matrix Convert a matrix solution from HDF5 to Matlab format.
+%
+% Signature:
+% hdf5_to_mat_matrix(infilename, solutionname, outfilename)
+%
+% Parameters:
+% infilename: the name of the HDF5 file containing the matrix solution
+% cellwise: true for a cell-wise solution vector, false for a vertex-wise
+% solutionname: the name of the solution vector inside the HDF5 file
+% outfilename: the name of the output Matlab .mat-file
+
+Points = h5read(infilename, '/Block_0_t000000/Geometry/Coordinates');
+Points = Points([1 2], :)';
+if cellwise
+  Pmixed = h5read(infilename, strcat('/Block_0_t000000/Cell/', solutionname));
+  P = [];
+end
+Tmixed = h5read(infilename, '/Block_0_t000000/Topology/DataArray0');
+mixedIdx = 1;
+eIdx = 1;
+T = [];
+while (mixedIdx < length(Tmixed))
+    eType = Tmixed(mixedIdx);
+    if eType == 4
+        T(eIdx, 1:3) = Tmixed(mixedIdx+1:mixedIdx+3)';
+        if cellwise
+            P(eIdx) = Pmixed(mixedIdx);
+        end
+        mixedIdx = mixedIdx + 4;
+        eIdx = eIdx + 1;
+    else
+        T(eIdx, 1:3) = Tmixed([mixedIdx+1, mixedIdx+2, mixedIdx+4])';
+        T(eIdx+1, 1:3) = Tmixed([mixedIdx+1, mixedIdx+4, mixedIdx+3])';
+        if cellwise
+            P(eIdx) = Pmixed(mixedIdx);
+            P(eIdx+1) = Pmixed(mixedIdx);
+        end
+        mixedIdx = mixedIdx + 5;
+        eIdx = eIdx + 2;
+    end
+end
+T = T + 1;
+if ~cellwise
+  pvertex = h5read(infilename, strcat('/Block_0_t000000/Node/', solutionname));
+  P = zeros(size(T, 1), 1);
+  factor = 1/size(T, 2);
+  for i = 1:size(T, 2)
+    P = P + factor.*pvertex(T(:, i));
+  end
+end
+
+save(outfilename, 'T', 'P', 'Points');
+fprintf('Wrote file %s.\n', outfilename);
+
+return

scripts/hdf5_to_mat_ref.m

+function hdf5_to_mat_ref(infilename, outfilename, domainarea)
+% hdf5_to_mat_ref Convert a reference solution from HDF5 to Matlab format.
+%
+% Signature:
+% hdf5_to_mat_ref(infilename, outfilename, domainarea)
+%
+% Parameters:
+% infilename: the name of the HDF5 file containing the matrix solution
+% outfilename: the name of the output Matlab .mat-file
+% domainarea: area of the computational domain
+
+tvec = h5read(infilename, '/Block_0_t000000/Topology');
+tvec = tvec' + 1;
+t = {};
+for i = 1:length(tvec)
+    t{i} = tvec(i, :);
+end
+t = t';
+x = h5read(infilename, '/Block_0_t000000/Geometry/Coordinates');
+x = double(x([1 2], :)');
+p = h5read(infilename, strcat('/Block_0_t000000/Cell/', 'pressure'));
+k = h5read(infilename, strcat('/Block_0_t000000/Cell/', 'permeability'));
+
+save(outfilename, 't', 'p', 'x', 'k', 'domainarea');
+fprintf('Wrote file %s.\n', outfilename);
+
+return

scripts/hdf5_to_mat_ref_mixed.m

+function hdf5_to_mat_ref_mixed(infilename, outfilename, domainarea)
+% hdf5_to_mat_ref_mixed Convert a reference solution from HDF5 to Matlab format.
+%
+% Signature:
+% hdf5_to_mat_ref_mixed(infilename, outfilename, domainarea)
+%
+% Parameters:
+% infilename: the name of the HDF5 file containing the matrix solution
+% outfilename: the name of the output Matlab .mat-file
+% domainarea: area of the computational domain
+
+x = h5read(infilename, '/Block_0_t000000/Geometry/Coordinates');
+x = x([1 2], :)';
+pmixed = h5read(infilename, strcat('/Block_0_t000000/Cell/', 'pressure'));
+p = [];
+kmixed = h5read(infilename, strcat('/Block_0_t000000/Cell/', 'permeability'));
+k = [];
+tmixed = h5read(infilename, '/Block_0_t000000/Topology/DataArray0');
+mixedIdx = 1;
+eIdx = 1;
+t = [];
+while (mixedIdx < length(tmixed))
+    eType = tmixed(mixedIdx);
+    if eType == 4
+        t(eIdx, 1:3) = tmixed(mixedIdx+1:mixedIdx+3)';
+        p(eIdx) = pmixed(mixedIdx);
+        k(eIdx) = kmixed(mixedIdx);
+        mixedIdx = mixedIdx + 4;
+        eIdx = eIdx + 1;
+    else
+        t(eIdx, 1:3) = tmixed([mixedIdx+1, mixedIdx+2, mixedIdx+4])';
+        t(eIdx+1, 1:3) = tmixed([mixedIdx+1, mixedIdx+4, mixedIdx+3])';
+        p(eIdx) = pmixed(mixedIdx);
+        p(eIdx+1) = pmixed(mixedIdx);
+        k(eIdx) = kmixed(mixedIdx);
+        k(eIdx+1) = kmixed(mixedIdx);
+        mixedIdx = mixedIdx + 5;
+        eIdx = eIdx + 2;
+    end
+end
+t = t + 1;
+
+save(outfilename, 't', 'p', 'x', 'k', 'domainarea');
+fprintf('Wrote file %s.\n', outfilename);
+
+return

scripts/line_norm.m

 function [relativeError, errorSquared, normalizationSquared] = ...
-                                    line_norm(fracXAndP, p, t, x, FRAC)
+                                    line_norm(fracXAndP, p, t, x, FRAC, maxp, minp)
 %line_norm calculates the error for a single line
 %   For the format of fracXAndP, see the Description file. p, t and x are
 %   the reference pressure, triangulation and point list, respectively.
@@ -92,7 +92,7 @@ end
 % norms in evaluate_norm:
 errorSquared = evaluate_norm(IntersectionPoints, intersectionPoints, POn, pOn, isVertical,endPoints);
 
-normalizationSquared = (max(p)-min(p))^2*fractureLength;
+normalizationSquared = (maxp-minp)^2*fractureLength;
 relativeError = sqrt((errorSquared)/(normalizationSquared));