From eeb491db76b4a99dcccfed13064f69fb85677df5 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Tue, 13 Jun 2017 23:56:06 +0200
Subject: [PATCH] [scripts] calculate the error only on the exact matrix domain
---
scripts/compute_1d_error.m | 9 ++++--
scripts/compute_2d_error.m | 7 +++--
scripts/hdf5_to_mat_mixed.m | 56 +++++++++++++++++++++++++++++++++
scripts/hdf5_to_mat_ref.m | 27 ++++++++++++++++
scripts/hdf5_to_mat_ref_mixed.m | 46 +++++++++++++++++++++++++++
scripts/line_norm.m | 4 +--
6 files changed, 143 insertions(+), 6 deletions(-)
create mode 100644 scripts/hdf5_to_mat_mixed.m
create mode 100644 scripts/hdf5_to_mat_ref.m
create mode 100644 scripts/hdf5_to_mat_ref_mixed.m
diff --git a/scripts/compute_1d_error.m b/scripts/compute_1d_error.m
index ec13cac..18a256e 100644
--- a/scripts/compute_1d_error.m
+++ b/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));
diff --git a/scripts/compute_2d_error.m b/scripts/compute_2d_error.m
index 4bab327..4a49a6a 100644
--- a/scripts/compute_2d_error.m
+++ b/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));
diff --git a/scripts/hdf5_to_mat_mixed.m b/scripts/hdf5_to_mat_mixed.m
new file mode 100644
index 0000000..f1e717a
--- /dev/null
+++ b/scripts/hdf5_to_mat_mixed.m
@@ -0,0 +1,56 @@
+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
diff --git a/scripts/hdf5_to_mat_ref.m b/scripts/hdf5_to_mat_ref.m
new file mode 100644
index 0000000..5cde848
--- /dev/null
+++ b/scripts/hdf5_to_mat_ref.m
@@ -0,0 +1,27 @@
+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
diff --git a/scripts/hdf5_to_mat_ref_mixed.m b/scripts/hdf5_to_mat_ref_mixed.m
new file mode 100644
index 0000000..283531d
--- /dev/null
+++ b/scripts/hdf5_to_mat_ref_mixed.m
@@ -0,0 +1,46 @@
+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
diff --git a/scripts/line_norm.m b/scripts/line_norm.m
index 8bbf0b5..bd6b624 100644
--- a/scripts/line_norm.m
+++ b/scripts/line_norm.m
@@ -1,5 +1,5 @@
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));
--
GitLab