From 84539c55432155f1c6ca46fd14c2d6f2e1fa7d78 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Thu, 15 Dec 2016 17:21:47 +0100
Subject: [PATCH] add Matlab scripts
---
scripts/EvaluateAllFractures.m | 32 ++++
scripts/Exact2dNorm.m | 117 ++++++++++++++
scripts/LineNorm.m | 127 +++++++++++++++
scripts/LineNormFucntion.m | 107 +++++++++++++
scripts/NormApproximation.m | 82 ++++++++++
scripts/area_on_coarse_point.m | 57 +++++++
scripts/area_on_face.m | 67 ++++++++
scripts/check_points.m | 18 +++
scripts/combine_faces.m | 23 +++
scripts/divideTriangulation.m | 59 +++++++
scripts/evaluate_norm.m | 192 +++++++++++++++++++++++
scripts/extract_on_line.m | 92 +++++++++++
scripts/intersection.m | 27 ++++
scripts/intersection_between_inclusive.m | 59 +++++++
scripts/intersections_of_cells.m | 40 +++++
15 files changed, 1099 insertions(+)
create mode 100644 scripts/EvaluateAllFractures.m
create mode 100644 scripts/Exact2dNorm.m
create mode 100644 scripts/LineNorm.m
create mode 100644 scripts/LineNormFucntion.m
create mode 100644 scripts/NormApproximation.m
create mode 100644 scripts/area_on_coarse_point.m
create mode 100644 scripts/area_on_face.m
create mode 100644 scripts/check_points.m
create mode 100644 scripts/combine_faces.m
create mode 100644 scripts/divideTriangulation.m
create mode 100644 scripts/evaluate_norm.m
create mode 100644 scripts/extract_on_line.m
create mode 100644 scripts/intersection.m
create mode 100644 scripts/intersection_between_inclusive.m
create mode 100644 scripts/intersections_of_cells.m
diff --git a/scripts/EvaluateAllFractures.m b/scripts/EvaluateAllFractures.m
new file mode 100644
index 0000000..c4d2386
--- /dev/null
+++ b/scripts/EvaluateAllFractures.m
@@ -0,0 +1,32 @@
+% Evaluates the norm in each of the fractures and adds the computed values
+% to the input file. This file must contain a cell array with fracxandp for
+% each of the fractures, see the Description for details on how each
+% fracxandp should be constructed.
+clear
+load geiger_permeable_reference.mat
+load anna_reference_lefttoright.mat
+% load geiger_blocking_reference.mat
+%load Hydrocoin_reference.mat
+%coarsefile = 'coarse file.mat';
+
+%coarsefile = '../hydrocoin/results/edfm_fractures.mat';
+coarsefile = '../geiger/results/permeable_edfm.mat';
+%coarsefile = '../geiger/results/blocking_edfm.mat';
+
+coarsefile = '/home/elle/Dropbox/2015-fractures-comparison/new_case_4/results_bis/edfm_fractures.mat';
+
+load(coarsefile)
+nfrac = length(fracxandp);
+Erel_f_vec = zeros(nfrac,1);
+E2_f = zeros(nfrac,1);
+dP2_f = zeros(nfrac,1);
+FRAC = true;
+%%
+for i = 1:nfrac
+ xandp = fracxandp{i};
+
+ [Erel_f_vec(i),E2_f(i),dP2_f(i)] = LineNormFucntion(xandp,p,t,x,FRAC);
+end
+Erel_f = sqrt(sum(E2_f)/sum(dP2_f))
+save(coarsefile,'Erel_f','E2_f','dP2_f','-append')
+
diff --git a/scripts/Exact2dNorm.m b/scripts/Exact2dNorm.m
new file mode 100644
index 0000000..5c89cc5
--- /dev/null
+++ b/scripts/Exact2dNorm.m
@@ -0,0 +1,117 @@
+%% Load
+clear
+tic
+coarse = '../new_case_4/results/bernd_dfm_lefttoright.mat';
+fine = 'anna_reference_lefttoright.mat';
+% coarse = 'permeable_kanskje_homo2.mat';
+% fine = 'geiger_permeable_reference.mat';
+load(coarse,'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(fine,'t','p','x','domainarea')
+
+points = x;
+
+
+nc = length(t);
+max_n_nodes = 4; % for the fine mesh
+t_nan = NaN(nc,max_n_nodes);
+remain = true(nc,1);
+facelength = zeros(nc,max_n_nodes);
+e = zeros(nc,1);
+plotE = zeros(nc,1);
+% on coarse corners
+for c = 1:length(t)
+ id = t{c};
+ t_nan(c,1:length(id)) = id;
+ pol = points(id,:);
+ [in,on] = inpolygon(Points(:,1),Points(:,2),pol(:,1),pol(:,2));
+ Point_id = find(in & ~on);
+ if ~isempty(Point_id)
+ if sum(Point_id)>1
+ uniquep = unique(1e-5.*round(1e5.*Points(Point_id,:)),'rows'); %some methods
+ % define certain points twice
+ if size(uniquep,1)>1
+ disp(uniquep);
+ error('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(Point_id)
+
+ Cs = find(any(T==Point_id(Pid),2));
+ for i = 1:length(Cs)
+ C = Cs(i);
+ Id = T(C,:);
+ Pol = Points(Id,:);
+ a = area_on_coarse_point(pol,Pol,Points(Point_id,:));
+ dp = P(C)-p(c);
+ e(c) = e(c) + (dp)^2*a;
+
+ afracdp = dp*a/polyarea(pol(:,1),pol(:,2));
+ plotE(c) = plotE(c) + afracdp;
+ end
+ end
+ remain(c) = false;
+
+ end
+end
+%%
+% eliminate the cells already taken care of
+t2 = t(remain);
+t_nan2 = t_nan(remain,:);
+nanind = isnan(t_nan2(:,4)); % obs bare for 3+4
+t3 = t_nan2;
+t3(nanind,4) = t_nan2(nanind,1);
+e2 = zeros(size(t_nan2,1),1);
+plotE2 = zeros(size(e2));
+p2 = p(remain);
+
+
+for C = 1:length(T)
+ % construct coarse cell
+ Id = T(C,:);
+ Pol = Points(Id,:);
+
+ % Find fine points inside the coarse cell
+ [in,on] = inpolygon(points(:,1),points(:,2),Pol(:,1),Pol(:,2));
+ points_in = find(in & ~on);
+
+ % Distinguish between cells completely and partly inside large cell
+ tr_in = ismember(t_nan2,points_in);
+ cells_in = find(any(tr_in,2));
+ cells_completely_in = cells_in(all(ismember(t3(cells_in,:),points_in),2));
+ cells_partly_in = setdiff(cells_in,cells_completely_in);
+
+ % Evaluate for the two cases
+ for i = 1:length(cells_completely_in)
+ c = cells_completely_in(i);
+ x = points(t2{c},:);
+ [~,a] = convhull(x);
+ dp = (P(C)-p2(c))^2;
+ e2(c) = dp*a;
+ plotE2(c) = P(C)-p2(c);
+ end
+ for i = 1:length(cells_partly_in)
+ c = cells_partly_in(i);
+ pol = points(t2{c},:);
+ a = area_on_face(pol,Pol);
+ dp = (P(C)-p2(c));
+ e2(c) = e2(c) + dp^2*a;
+ afracdp = dp*a/polyarea(pol(:,1),pol(:,2));
+ plotE2(c) = plotE2(c) + afracdp;
+ end
+end
+e(remain) = e2;
+plotE(remain) = plotE2;
+
+toc
+deltaP = max(p)-min(p);
+Erel_m = sqrt(sum(e))/(deltaP*sqrt(domainarea));
+save(coarse,'Erel_m','plotE','-append')
diff --git a/scripts/LineNorm.m b/scripts/LineNorm.m
new file mode 100644
index 0000000..414aa91
--- /dev/null
+++ b/scripts/LineNorm.m
@@ -0,0 +1,127 @@
+coarsefile = '../new_case_4/results/bernd_dfm_lefttoright.mat';
+load(coarsefile);
+%load 'coarse file.mat'
+X = Points;
+
+FRAC = true;
+if FRAC
+% i = 2; %fracture identity
+ Xp = fracxandp{i}; % The fracture information has here the format
+ endpoints=[Xp(1, 1:2); Xp(end, 1:2)];
+
+ % [x1(c) y1(c) p(c); x2(c) y2(c) p(c)] for cell c in position
+ % Xp(2*c-1:2*c, :), the coordinates being the vertices of the fracture.
+end
+
+%load geiger_permeable_reference.mat
+load anna_reference_lefttoright.mat
+
+
+% Both the "coarse" and "fine" input files are assumed to have a triangulation
+% CELL t with one (number of vertices of the cell) x 1 pointer to the points
+% x of each grid cell and vector p for pressure in each cell (length(p) ==
+% length(t)). In case a fracture is being evaluated, the only coarse
+% simulation information needed is fracxandp, a cell containing the
+% pressure and coordinate information for each fracture.
+% Capital Letters Indicate The Coarse Grid.
+
+boundary_line = false;
+
+%segment = [endpoints(1),endpoints(2,1)];
+segment = [1, 1];
+% assign two identical values for the whole cross-section of
+% the domain. Assign x values for the endpoints of the desired segment
+% unless the line is vertical, in which case the two y values should be
+% provided.
+% NOTE the segments are assumed to be chosen so that the cells end at the
+% endpoints (fracture endpoints, in our case). If not, only the part of the
+% line that is covered by coarse and fine cells lying in the interior of
+% the segment is evaluated.
+if endpoints(1)==endpoints(2)
+ a=endpoints(1);
+ b=0;
+ isvertical = true;
+ if segment(1) == segment(2)
+ fracture_length = abs(endpoints(3)-endpoints(4));
+ else
+ fracture_length = abs(segment(1)-segment(2));
+ end
+else
+ a = (endpoints(2,2) - endpoints(1,2)) / (endpoints(2,1)-endpoints(1,1));
+ b = endpoints(1,2) - a*endpoints(1,1);
+ isvertical = false;
+ if segment(1) == segment(2)
+ fracture_length = hypot(abs(endpoints(1)-endpoints(2)), ...
+ abs(endpoints(3)-endpoints(4)));
+ else
+ fracture_length = hypot(abs(segment(1)-segment(2)), ...
+ abs(endpoints(3)-endpoints(4)));
+ end
+end
+
+
+[nw, vertices_on] = check_points(x,a,b,isvertical);
+if ~FRAC
+ [Nw, Vertices_on] = check_points(X,a,b,isvertical);
+end
+%must unfortunately allow for cells with different number of vertices both
+%for fine grid (Bernd uses both for the hydrocoin case, at least) and coarse.
+% Therefore, cell_points_on is a CELL with nvertices(c) pointer extracted
+%from t. c_p_nw_on is the corresponding logical nw cell.
+
+[cell_points_on, p_on, c_p_nw_on, cell_vertices_on, vertices_p] = ...
+ extract_on_line(t,x,nw,vertices_on,p);
+if ~FRAC
+ [Cell_points_on, P_on, C_p_nw_on, Cell_vertices_on, Vertices_p] = ...
+ extract_on_line(T,X,Nw,Vertices_on,P);
+else
+ Cell_vertices_on = [];
+ Vertices_p = [];
+end
+
+
+if ~boundary_line && ~isempty(cell_vertices_on)
+
+ [cell_vertices_on, vertices_p] = combine_faces(cell_vertices_on, vertices_p,isvertical);
+end
+if ~boundary_line && ~isempty(Cell_vertices_on)
+
+ [Cell_vertices_on, Vertices_p] = combine_faces(Cell_vertices_on, Vertices_p,isvertical);
+end
+
+
+% find the points where the line intersects each of the fine and coarse
+% cells, respectively. Format intersectionpoints(c) = [x1,y1,x2,y2] for the
+% two points for cell c.
+
+
+[intersectionpoints] = intersections_of_cells(endpoints, isvertical, ...
+ cell_points_on, c_p_nw_on);
+
+if FRAC
+ n = length(Xp);
+ ind_two = linspace(2,n,n/2);
+ ind_one = ind_two-1;
+ Intersectionpoints = [Xp(ind_one,1:2), Xp(ind_two,1:2)];
+ P_on = Xp(ind_two,3);
+else
+ [Intersectionpoints] = intersections_of_cells(endpoints, isvertical, ...
+ Cell_points_on, C_p_nw_on);
+end
+
+% Add the values for the faces coinciding with the line:
+if ~isempty(cell_vertices_on)
+ intersectionpoints = [intersectionpoints;cell_vertices_on];
+ p_on = [p_on;vertices_p];
+end
+if ~isempty(Cell_vertices_on)
+ Intersectionpoints = [Intersectionpoints;Cell_vertices_on];
+ P_on = [P_on;Vertices_p];
+end
+
+
+% Loop through large cells, find smalls cells (partly) inside and evaluate
+% norms
+[E2] = evaluate_norm(Intersectionpoints, intersectionpoints, P_on, p_on, isvertical,segment);
+dP2 = (max(p)-min(p))^2*fracture_length;
+Erel = sqrt(E2)/sqrt(dP2)
diff --git a/scripts/LineNormFucntion.m b/scripts/LineNormFucntion.m
new file mode 100644
index 0000000..2c0fecf
--- /dev/null
+++ b/scripts/LineNormFucntion.m
@@ -0,0 +1,107 @@
+function [Erel,E2,dP2] = LineNormFucntion(Xp,p,t,x,FRAC)
+%UNTITLED2 Summary of this function goes here
+% Detailed explanation goes here
+
+
+
+
+boundary_line = false;
+
+% assign two identical values for the whole cross-section of
+% the domain. Assign x values for the endpoints of the desired segment
+% unless the line is vertical, in which case the two y values should be
+% provided.
+% NOTE the segments are assumed to be chosen so that the cells end at the
+% endpoints (fracture endpoints, in our case). If not, only the part of the
+% line that is covered by coarse and fine cells lying in the interior of
+% the segment is evaluated.
+if Xp(1)==Xp(2)
+ endpoints = [Xp(1),min(Xp(:,2)); ...
+ Xp(1),max(Xp(:,2))];
+ a=endpoints(1);
+ b=0;
+ isvertical = true;
+ segment = endpoints(3:4);
+ fracture_length = endpoints(4)-endpoints(3);
+else
+ [x1,i_min] = min(Xp(:,1));
+ [x2,i_max] = max(Xp(:,1));
+ endpoints = [x1,Xp(i_min,2);x2,Xp(i_max,2)];
+ a = (endpoints(2,2) - endpoints(1,2)) / (endpoints(2,1)-endpoints(1,1));
+ b = endpoints(1,2) - a*endpoints(1,1);
+ isvertical = false;
+ segment = endpoints(1:2);
+ fracture_length = hypot(abs(endpoints(1)-endpoints(2)), abs(endpoints(3)-endpoints(4)));
+end
+
+
+[nw, vertices_on] = check_points(x,a,b,isvertical);
+if ~FRAC
+ [Nw, Vertices_on] = check_points(X,a,b,isvertical);
+end
+%must unfortunately allow for cells with different number of vertices both
+%for fine grid (Bernd uses both for the hydrocoin case, at least) and coarse.
+% Therefore, cell_points_on is a CELL with nvertices(c) pointer extracted
+%from t. c_p_nw_on is the corresponding logical nw cell.
+
+[cell_points_on, p_on, c_p_nw_on, cell_vertices_on, vertices_p] = ...
+ extract_on_line(t,x,nw,vertices_on,p);
+if ~FRAC
+ [Cell_points_on, P_on, C_p_nw_on, Cell_vertices_on, Vertices_p] = ...
+ extract_on_line(T,X,Nw,Vertices_on,P);
+else
+ Cell_vertices_on = [];
+ Vertices_p = [];
+end
+
+
+if ~boundary_line && ~isempty(cell_vertices_on)
+
+ [cell_vertices_on, vertices_p] = combine_faces(cell_vertices_on, vertices_p,isvertical);
+end
+if ~boundary_line && ~isempty(Cell_vertices_on)
+
+ [Cell_vertices_on, Vertices_p] = combine_faces(Cell_vertices_on, Vertices_p,isvertical);
+end
+
+
+% find the points where the line intersects each of the fine and coarse
+% cells, respectively. Format intersectionpoints(c) = [x1,y1,x2,y2] for the
+% two points for cell c.
+
+
+[intersectionpoints] = intersections_of_cells(endpoints, isvertical, ...
+ cell_points_on, c_p_nw_on);
+
+if FRAC
+ n = size(Xp,1);
+ ind_two = linspace(2,n,n/2);
+ ind_one = ind_two-1;
+ Intersectionpoints = [Xp(ind_one,1:2), Xp(ind_two,1:2)];
+ P_on = Xp(ind_two,3);
+else
+ [Intersectionpoints] = intersections_of_cells(endpoints, isvertical, ...
+ Cell_points_on, C_p_nw_on);
+end
+
+% Add the values for the faces coinciding with the line:
+if ~isempty(cell_vertices_on)
+ intersectionpoints = [intersectionpoints;cell_vertices_on];
+ p_on = [p_on;vertices_p];
+end
+if ~isempty(Cell_vertices_on)
+ Intersectionpoints = [Intersectionpoints;Cell_vertices_on];
+ P_on = [P_on;Vertices_p];
+end
+
+
+% Loop through large cells, find smalls cells (partly) inside and evaluate
+% norms
+
+E2 = evaluate_norm(Intersectionpoints, intersectionpoints, P_on, p_on, isvertical,segment);
+dP2 = (max(p)-min(p))^2*fracture_length;
+Erel = sqrt((E2)/(dP2));
+
+
+end
+
diff --git a/scripts/NormApproximation.m b/scripts/NormApproximation.m
new file mode 100644
index 0000000..b086de4
--- /dev/null
+++ b/scripts/NormApproximation.m
@@ -0,0 +1,82 @@
+clear;
+%coarsefile = '../hydrocoin/results/edfm.mat';
+coarsefile = '../geiger/results/permeable_edfm.mat';
+%coarsefile = '../geiger/results/blocking_edfm.mat';
+
+load(coarsefile);
+% Should contain a pressure vector P of length equal to the number of cells
+% with dimension two (lower-dimensional fracture cells need not be included,
+% but will not affect the result).
+% If the triangulation uses cells with different number of vertices, it should
+% be a cell array T. If not, a matrix T with size(T) = [ # cells,
+% numberofvertices/cell] works just as well. The triangulation points
+% should be provided in the matrix Points, size(Points) = [ # points, 2 ].
+
+% load geiger_blocking_reference.mat
+%load Hydrocoin_reference.mat
+load geiger_permeable_reference.mat
+
+% Contains a cell center, area and pressure vector, each with length equal
+% to the number of cells in the reference solution.
+
+
+
+% in the following, variables starting with Capital Letters Are "Coarse"!
+
+
+X = Points(:,1);
+Y = Points(:,2);
+
+
+%Ea is a vector and count is introduced to help analyze the error evaluation
+% (e.g. max(Ea) >> max(Ea\max(Ea)) could be suspicious, and count ~= ncells
+% most certainly is). If the runtime is too high, Ea could be replaced by
+% a scalar and count left out.
+Ea = zeros(size(p));
+% normalization = zeros(size(p));
+count=0;
+
+% Outer loop for the coarse cells:
+for C = 1:length(P)
+ %in contains the fine cells with centre inside the coarse cell C
+ if iscell(T)
+ Cellpoints = T{C};
+ else
+ Cellpoints = T(C,:);
+ end
+
+ X_vec = X(Cellpoints); Y_vec = Y(Cellpoints);
+ hull = convhull( X_vec, Y_vec );
+ hull = hull( 1:(end-1) );
+
+ in = inpolygon(cc(:,1), cc(:,2), X_vec( hull ), Y_vec( hull ) );
+ in = find(in);
+
+ %Evaluate the error on the fine grid:
+ for i = 1:length(in)
+ c = in(i);
+ Ea(c) = Ea(c)+ a(c)*(p(c)-P(C))^2;
+% normalization(c) = normalization(c) + a(c)*p(c)^2;
+ count = count+1;
+ end
+end
+
+normalization = a*(max(p)-min(p))^2;
+E2_m = sum(Ea);
+dP2_m = sum(normalization);
+Erel_m = sqrt(E2_m/dP2_m)
+
+
+% To save some time (relevant on my computer, at least), only a selection
+% of points is plotted. If one wishes to see all points, simply change
+% nplotpoints to length(Ea). Choose between the scatter3 version you find
+% most revealing.
+nplotpoints = 3e4;
+ind = floor(linspace(1, length(normalization), nplotpoints));
+eloc = sqrt(Ea(ind))./sqrt(normalization(ind));
+figure
+ scatter3(cc(ind,1),cc(ind,2),eloc);
+%scatter3(cc(ind,1),cc(ind,2),eloc,1,eloc)
+%scatter3(cc(ind,1),cc(ind,2),zeros(size(ind)),1,eloc)
+
+save(coarsefile,'Erel_m','E2_m','dP2_m','-append')
diff --git a/scripts/area_on_coarse_point.m b/scripts/area_on_coarse_point.m
new file mode 100644
index 0000000..951c944
--- /dev/null
+++ b/scripts/area_on_coarse_point.m
@@ -0,0 +1,57 @@
+function [ a ] = area_on_coarse_point( pol,Pol,Point )
+% Calculates intersection areas of convex polygons
+% For a small polygon surrounding a point in a much coarser triangulation,
+% this function calculates the intersection area with the coarse polygon
+% Pol. Point is one of the points defining Pol. Note that both size
+% difference and convexity assumptions must hold!
+% f1 er point 1 og point 2, f2 er p2 og p3 ... fn er pn og p1
+N = size(Pol,1);
+n = size(pol,1);
+inPol = find(ismember(Pol,Point,'rows'));
+
+% Find the faces of Pol adjacent to Point
+if inPol == 1
+ Firstface = [Pol(N,:);Pol(1,:)];
+else
+ Firstface = [Pol(inPol-1,:);Pol(inPol,:)];
+end
+if inPol == N
+ Secondface = [Pol(N,:);Pol(1,:)];
+else
+ Secondface = [Pol(inPol,:);Pol(inPol+1,:)];
+end
+
+pol2 = [pol;pol(1,:)];
+for i = 1:n
+ face = [pol2(i,:);pol2(i+1,:)];
+ intersection1 = intersection_between_inclusive(Firstface, ...
+ Firstface(1)==Firstface(2), face, face(1)==face(2),true);
+ if ~isempty(intersection1)
+ break
+ end
+end
+for j = 1:n
+ face = [pol2(j,:);pol2(j+1,:)];
+ intersection2 = intersection_between_inclusive(Secondface, ...
+ Secondface(1)==Secondface(2), face, face(1)==face(2),true);
+ if ~isempty(intersection2)
+ break
+ end
+end
+if i ~= j
+ smallpointsinside = inpolygon(pol(:,1),pol(:,2),Pol(:,1),Pol(:,2));
+ smallpointsinside = pol(find(smallpointsinside),:);
+else
+ smallpointsinside = [];
+end
+newpoints = [Point; intersection1;intersection2;smallpointsinside];
+
+
+if size( unique( [ newpoints(:,1) newpoints(:,2) ], 'rows' ), 1 ) > 2
+ [~, a] = convhull(newpoints(:,1),newpoints(:,2));
+else
+ a = 0;
+end
+
+end
+
diff --git a/scripts/area_on_face.m b/scripts/area_on_face.m
new file mode 100644
index 0000000..fa721b1
--- /dev/null
+++ b/scripts/area_on_face.m
@@ -0,0 +1,67 @@
+function [ a ] = area_on_face( pol,Pol )
+%area_on_face calculates intersection area between convex polygons
+% The small polygon pol is cut by one of the edges of the much larger
+% Pol. Their intersection area is computed.
+N = size(Pol,1);
+n = size(pol,1);
+nw = false(n,N);
+on = false(n,N);
+alloneside = false(N,1);
+Pol2 = [Pol;Pol(1,:)];
+pol2 = [pol;pol(1,:)];
+isv = false(N,1);
+for i = 1:N
+ if Pol2(i)==Pol2(i+1)
+ isv(i) = true;
+ a = Pol2(i);
+ b = 0;
+ else
+ a = (Pol2(i+1,2)-Pol2(i,2))/(Pol2(i+1)-Pol2(i));
+ b = Pol2(i,2)-a*Pol2(i);
+ end
+ [nw(:,i),on(:,i)] = check_points(pol,a,b,isv(i));
+ alloneside(i) = all(nw(:,i)) || all(~nw(:,i));
+end
+
+% [nw(:,N),on(:,N)] = check_points(pol,a,b,isv); % hva med på linje??
+% alloneside(N) = all(nw(:,N)) || all(~nw(:,N));
+
+
+Faces = find(~alloneside);
+Nf = length(Faces);
+ic = zeros(Nf*2,2);
+j = 0;
+nw2 = [nw;nw(1,:)];
+on2 = [on;on(1,:)];
+for k = 1:Nf
+ F = Faces(k);
+ for i = 1:n
+ if on(i,F)
+ if on2(i+1,F)
+ s = sortrows([pol2(i:i+1,:);Pol2(F:F+1,:)]);
+ ic(j+1:j+2,:) = s(2:3,:);
+ j = j+2;
+ else
+ j = j +1;
+ ic(j,:) = pol(i,:);
+ end
+ elseif on2(i+1,F)
+
+ elseif nw2(i,F) ~=nw2(i+1,F)
+ j = j+1;
+ ic2 = intersection_between_inclusive([Pol2(F),Pol2(F,2);Pol2(F+1),Pol2(F+1,2)], ...
+ isv(F),[pol2(i),pol2(i,2);pol2(i+1),pol2(i+1,2)],pol2(i)==pol2(i+1),true);
+ if ~isempty(ic2)
+ ic(j,:) = ic2;
+ else
+ j = j - 1;
+ end
+ end
+ end
+end
+
+
+inside = inpolygon(pol(:,1),pol(:,2),Pol(:,1),Pol(:,2));
+insidepoints = pol(inside,:);
+[~, a] = convhull([ic(1:j,:);insidepoints]);
+end
\ No newline at end of file
diff --git a/scripts/check_points.m b/scripts/check_points.m
new file mode 100644
index 0000000..7e9bb65
--- /dev/null
+++ b/scripts/check_points.m
@@ -0,0 +1,18 @@
+function [ nw, on] = check_points( x, a,b , isvertical)
+%check_points Determines which of a number of 2D points in x lie north (west)
+% of the line ax + b and if any points lie on the line
+% If the line is vertical, the constant x-value goes in the input a
+
+if isvertical
+ nw = x(:,1) < a;
+ on = x(:,1) == a;
+elseif a == 0
+ nw = x(:,2) > b;
+ on = x(:,2) == b;
+else
+ nw = a*x(:,1) + b < x(:,2);
+ on = a*x(:,1) + b == x(:,2);
+end
+
+end
+
diff --git a/scripts/combine_faces.m b/scripts/combine_faces.m
new file mode 100644
index 0000000..c77bd9d
--- /dev/null
+++ b/scripts/combine_faces.m
@@ -0,0 +1,23 @@
+function [cell_vertices_on, vertices_p] = combine_faces(cell_vertices_on, vertices_p,isvertical)
+%UNTITLED2 calculates the average of the two pressure values for the faces
+% which coincide with the line
+
+if isvertical
+ cell_xory = sort([cell_vertices_on(:,2),cell_vertices_on(:,4),],2);
+else
+ cell_xory = sort([cell_vertices_on(:,1),cell_vertices_on(:,3),],2);
+end
+for i = 1:length(vertices_p)/2 % two values for each face
+ ind = ismember(cell_xory(:,1),cell_xory(i,1));
+ r = find(ind);
+ n = length(ind);
+ % remove the second cell
+ cell_vertices_on = [cell_vertices_on(1:r(2)-1,:);cell_vertices_on(r(2)+1:n,:)];
+ cell_xory = [cell_xory(1:r(2)-1,:); cell_xory(r(2)+1:n, :)];
+ % replace the first pressure by the mean of the two
+ vertices_p(r(1)) = mean(vertices_p(ind));
+ % and remove the second one
+ vertices_p = [vertices_p(1:r(2)-1);vertices_p(r(2)+1:n)];
+
+end
+end
diff --git a/scripts/divideTriangulation.m b/scripts/divideTriangulation.m
new file mode 100644
index 0000000..94384fe
--- /dev/null
+++ b/scripts/divideTriangulation.m
@@ -0,0 +1,59 @@
+load SNPI_Bernd_raw.mat;
+% saved in importerref.m. Contains points, pressure and triangulation,
+% named x, p and t respectively.
+
+f = t(:,4)~=0;
+th = ~f;
+
+% last letter "t" for three, "f" for four
+pt = p(th);
+pf = p(f);
+
+
+toldt = t(th, 1:3);
+toldf = t(f, :);
+[uniquet, intoldt, inuniquet] = unique(toldt);
+[uniquef, ~, inuniquef] = unique(toldf);
+
+xt = x(uniquet,:);
+xf = x(uniquef,:);
+
+% tnewt = zeros(size(toldt));
+% tnewf = zeros(size(toldf));
+%these sets of loops might perhaps be avoided by some "unique-trick". Might
+%reduce run time significantly
+% for i = 1:length(toldt)
+% for j = 1:3
+% tnewt(i,j)=find(uniquet==toldt(i,j));
+% end
+% end
+tnewt = reshape(inuniquet, size(toldt));
+tnewf = reshape(inuniquef, size(toldf));
+% for i = 1:length(toldf)
+% for j = 1:4
+% tnewf(i,j)=find(uniquef==toldf(i,j));
+% end
+% end
+%tt = triangulation(tnewt,xt);
+%can't use triangulation on the quads, so connectivity list and points are
+% saved seperately.
+clF = tnewf;
+%save('E:\Comparison paper\Mortar\SNPI_Bref', ...
+ % 'clF', 'xf', 'tt', 'pf', 'pt');
+
+
+
+Xt = xt(tnewt);
+yt = xt(:,2);
+Yt = yt(tnewt);
+Xf = xf(tnewf);
+yf = xf(:,2);
+Yf = yf(tnewf);
+midt = [sum(Xt,2)/3, sum(Yt,2)/3];
+areat = polyarea(Xt, Yt, 2);
+midf = [sum(Xf,2)/4, sum(Yf,2)/4];
+areaf = polyarea(Xf, Yf, 2);
+cc = [midt;midf];
+a = [areat;areaf];
+p = [pt; pf];
+save('C:\Users\Ivar\Desktop\L2norm\SNPI_reference_values', 'cc','a','p')
\ No newline at end of file
diff --git a/scripts/evaluate_norm.m b/scripts/evaluate_norm.m
new file mode 100644
index 0000000..15be47b
--- /dev/null
+++ b/scripts/evaluate_norm.m
@@ -0,0 +1,192 @@
+function [ E_abs2 ] = evaluate_norm( Intersectionpoints, intersectionpoints, ...
+ P, p, isvertical,segment)
+%UNTITLED Summary of this function goes here
+% Returns L2norm of error and fine pressure, both squared.
+ishorizontal = abs(Intersectionpoints(1,2) - Intersectionpoints(1,4))<eps;
+isp = intersectionpoints;
+Isp = Intersectionpoints;
+%sorting the points to have smallest point first in each
+%intersectionpointsline
+if isvertical
+ for c = 1:length(p)
+ if isp(c,2)>isp(c,4)
+ intersectionpoints(c,:) = [isp(c,3:4),isp(c,1:2)];
+ end
+ end
+ for C = 1:length(P)
+ if Isp(C,2)>Isp(C,4)
+ Intersectionpoints(C,:) = [Isp(C,3:4),Isp(C,1:2)];
+ end
+ end
+
+else
+ for c = 1:length(p)
+ if isp(c,1)>isp(c,3)
+ intersectionpoints(c,:) = [isp(c,3:4),isp(c,1:2)];
+ end
+ end
+ for C = 1:length(P)
+ if Isp(C,1)>Isp(C,3)
+ Intersectionpoints(C,:) = [Isp(C,3:4),Isp(C,1:2)];
+ end
+ end
+
+end
+ascending = sortrows([intersectionpoints,p]);
+Ascending = sortrows([Intersectionpoints,P]);
+if segment(1) ~= segment(2)
+ if isvertical
+ first = find(ascending(:,2) >= segment(1),1);
+ last = find(ascending(:,4) <= segment(2),1,'last');
+ ascending = ascending(first:last,:);
+
+ First = find(Ascending(:,2) >= segment(1),1);
+ Last = find(Ascending(:,4) <= segment(2),1,'last');
+ Ascending = Ascending(First:Last,:);
+ else
+ first = find(ascending(:,1) >= segment(1),1);
+ last = find(ascending(:,3) <= segment(2),1,'last');
+ ascending = ascending(first:last,:);
+
+ First = find(Ascending(:,1) >= segment(1),1);
+ Last = find(Ascending(:,3) <= segment(2),1,'last');
+ Ascending = Ascending(First:Last,:);
+
+ end
+end
+intersectionpoints = ascending(:,1:4);
+p = ascending(:,5);
+Intersectionpoints = Ascending(:,1:4);
+P = Ascending(:,5);
+
+E_abs2 = 0;
+
+for C = 1:length(P)
+ if isvertical
+ %get the [lower, upper] coarse points
+ Yc = sort ([Intersectionpoints(C,2),Intersectionpoints(C,4)]);
+
+
+ larger = intersectionpoints(:,4)>Yc(1);
+ smaller = intersectionpoints(:,2)<Yc(2);
+ are_in = find(larger + smaller == 2); % all the relevant fine cells,
+ % pointing to location in intersectionpoints (and local "p" =
+ % global "p_on")
+
+ % only one point of the large cell on the line
+ if isempty(are_in)
+
+
+ % fine cell is larger (on the line) than coarse cell
+ elseif length(are_in) == 1
+ larger_length = Yc(2)-Yc(1);
+ E_abs2 = E_abs2 + larger_length*(p(are_in)-P(C))^2;
+
+ % at least two fine cells (partly) inside the coarse
+ else
+ bottom = are_in(1);
+ top = are_in(end);
+ bottom_length = intersectionpoints(bottom,4) - Yc(1);
+ top_length = Yc(2) - intersectionpoints(top,2);
+ E_abs2 = E_abs2 + sum([bottom_length,top_length].* ...
+ ([p(bottom),p(top)]-P(C)).^2);
+
+
+ completely_in = are_in(2:end-1);
+ if ~isempty(completely_in)
+ completely_length = intersectionpoints(completely_in, 4) ...
+ -intersectionpoints(completely_in,2);
+ E_abs2 = E_abs2 + sum(completely_length.*(p(completely_in)-P(C)).^2);
+ end
+ end
+
+
+ %same cheap length evaluation applys for horizontal case
+ elseif ishorizontal
+ %get the [lower, upper] coarse points
+ Xc = sort ([Intersectionpoints(C,1),Intersectionpoints(C,3)]);
+
+
+ larger = intersectionpoints(:,3)>Xc(1);
+ smaller = intersectionpoints(:,1)<Xc(2);
+ are_in = find(larger + smaller == 2); % all the relevant fine cells,
+ % pointing to location in intersectionpoints (and local "p" =
+ % global "p_on")
+
+ % only one point of the large cell on the line
+ if isempty(are_in)
+
+
+ % fine cell is larger (on the line) than coarse cell
+ elseif length(are_in) == 1
+ larger_length = Xc(2)-Xc(1);
+ E_abs2 = E_abs2 + larger_length*(p(are_in)-P(C))^2;
+ % at least two fine cells (partly) inside the coarse
+ else
+ bottom = are_in(1);
+ top = are_in(end);
+ bottom_length = intersectionpoints(bottom,3) - Xc(1);
+ top_length = Xc(2) - intersectionpoints(top,1) ;
+ E_abs2 = E_abs2 + sum([bottom_length,top_length].* ...
+ ([p(bottom),p(top)]-P(C)).^2);
+
+
+ completely_in = are_in(2:end-1);
+ if ~isempty(completely_in)
+ completely_length = intersectionpoints(completely_in, 3) ...
+ -intersectionpoints(completely_in,1);
+ E_abs2 = E_abs2 + sum(completely_length.*(p(completely_in)-P(C)).^2);
+ end
+ end
+
+
+
+
+
+ % general, non-constant line. The following code also works on horizontal,
+ % but not vertical lines. The computation is a bit more expensive,
+ % hence the "elseif ishorizontal" above
+ else
+ %get the [lower, upper] coarse points
+ Xc2 = sortrows([Intersectionpoints(C,1:2);Intersectionpoints(C,3:4)]);
+ Xc = [Xc2(1:2)];
+ Yc = [Xc2(3:4)];
+
+ larger = intersectionpoints(:,3)>Xc(1);
+ smaller = intersectionpoints(:,1)<Xc(2);
+ are_in = find(larger + smaller == 2); % all the relevant fine cells,
+ % pointing to location in intersectionpoints (and local "p" =
+ % global "p_on")
+
+ % only one point of the large cell on the line
+ if isempty(are_in)
+
+
+ % fine cell is larger (on the line) than coarse cell
+ elseif length(are_in) == 1
+ larger_length = norm(Intersectionpoints(C,1:2) - ...
+ Intersectionpoints(C,3:4), 2);
+ E_abs2 = E_abs2 + larger_length*(p(are_in)-P(C))^2;
+ % at least two fine cells (partly) inside the coarse
+ else
+ bottom = are_in(1);
+ top = are_in(end);
+
+ bottom_length = norm(intersectionpoints(bottom,3:4) - ...
+ [Xc(1), Yc(1)], 2);
+ top_length = norm(intersectionpoints(top,1:2) - ...
+ [Xc(2), Yc(2)], 2);
+ E_abs2 = E_abs2 + sum([bottom_length,top_length].* ...
+ ([p(bottom),p(top)]-P(C)).^2);
+
+ completely_in = are_in(2:end-1);
+ if ~isempty(completely_in)
+ completely_length = hypot(intersectionpoints(completely_in, 1) ...
+ -intersectionpoints(completely_in,3), ...
+ intersectionpoints(completely_in, 2) - ...
+ intersectionpoints(completely_in,4));
+ E_abs2 = E_abs2 + sum(completely_length.*(p(completely_in)-P(C)).^2);
+ end
+ end
+ end
+end
\ No newline at end of file
diff --git a/scripts/extract_on_line.m b/scripts/extract_on_line.m
new file mode 100644
index 0000000..55b0d07
--- /dev/null
+++ b/scripts/extract_on_line.m
@@ -0,0 +1,92 @@
+function [ cell_points_on, p_on, cell_points_nw_on, cell_face_on,p_on2 ] ...
+ = extract_on_line( t, x, nw,on, p )
+%UNTITLED4 returns cell_points (pointers to x) and p for each cell
+%on the line. cell_on_line is the logical vector s.t. p(c_o_l)=p_on.
+%cell_points_nw_on has a cell for each grid cell with the logical information
+%whether the cell_points_nw_on{c}(i) pointed to by cell_points{c}(i) lies nw
+% NB for generality, t is assumed to be a cell. Note that in the case
+% that a face lies on the line, a matrix is filled, because we know the
+% exact number of points (2) in this case.
+
+count=1;
+count2=1;
+
+cell_points_on = cell(0,1);
+cell_points_nw_on = cell(0,1);
+cell_face_on = zeros(0,4);
+p_on = zeros(0,1);
+p_on2 = zeros(0,1);
+noton = ~on;
+if iscell(t);
+
+ for c = 1:length(t)
+ cells = t{c};
+ points_nw = nw(cells);
+ cell_on_line = not(all(points_nw) | all(not(points_nw)));
+ vert_on = on(cells);
+ % two vertices on the line
+ if sum(vert_on) == 2
+ vert_not = noton(cells);
+ not_on_nw = points_nw(vert_not);
+ % check that it is really a face, not e.g. point 1 and 3 of a
+ % quadrilateral
+ if all(not_on_nw) || all(~not_on_nw)
+ ind = find(vert_on);
+ cell_face_on(count2,1:2) = x(cells(ind(1)),:);
+ cell_face_on(count2,3:4) = x(cells(ind(2)),:);
+ p_on2(count2,1) = p(c);
+ count2 = count2+1;
+ %else it should be included in the main cell:
+ else
+ cell_points_on{count} = x(cells,:);
+ cell_points_nw_on{count} = points_nw;
+ p_on(count,1) = p(c);
+ count = count + 1;
+ end
+ % cells intersected by the line, not two vertices on it
+ elseif cell_on_line
+ cell_points_on{count} = x(cells,:);
+ cell_points_nw_on{count} = points_nw;
+ p_on(count,1) = p(c);
+ count = count + 1;
+
+ end
+ end
+else
+ for c = 1:length(t)
+ cells = t(c,:);
+ points_nw = nw(cells);
+ cell_on_line = not(all(points_nw) | all(not(points_nw)));
+ vert_on = on(cells);
+ % two vertices on the line
+ if sum(vert_on) == 2
+ vert_not = noton(cells);
+ not_on_nw = points_nw(vert_not);
+ % check that it is really a face, not e.g. point 1 and 3 of a
+ % quadrilateral
+ if all(not_on_nw) || all(~not_on_nw)
+ ind = find(vert_on);
+ cell_face_on(count2,1:2) = x(cells(ind(1)),:);
+ cell_face_on(count2,3:4) = x(cells(ind(2)),:);
+ p_on2(count2,1) = p(c);
+ count2 = count2+1;
+ %else it should be included in the main cell:
+ else
+ cell_points_on{count} = x(cells,:);
+ cell_points_nw_on{count} = points_nw;
+ p_on(count,1) = p(c);
+ count = count + 1;
+ end
+ % cells intersected by the line, not two vertices on it
+ elseif cell_on_line
+ cell_points_on{count} = x(cells,:);
+ cell_points_nw_on{count} = points_nw;
+ p_on(count,1) = p(c);
+ count = count + 1;
+
+ end
+ end
+end
+
+end
+
diff --git a/scripts/intersection.m b/scripts/intersection.m
new file mode 100644
index 0000000..d083df9
--- /dev/null
+++ b/scripts/intersection.m
@@ -0,0 +1,27 @@
+function [ intersection ] = intersection( k, vk, l, vl)
+%UNTITLED2 Calculates intersections of two lines specified by two sets of
+%points
+% lines in format [x1, y1; x2, y2]
+% [11, 12; 21, 22]
+% line k = ax + b, l = cx + d
+% vi = true iff i vertical
+if vk
+ c = (l(2,2) - l(1,2)) / (l(2,1)-l(1,1));
+ d = l(1,2) - c*l(1,1);
+ intersection = [k(1), c * k(1) + d];
+elseif vl
+ a = (k(2,2) - k(1,2)) / (k(2,1)-k(1,1));
+ b = k(1,2) - a*k(1,1);
+ intersection = [l(1), a * l(1) + b];
+else
+ a = (k(2,2) - k(1,2)) / (k(2,1)-k(1,1));
+ c = (l(2,2) - l(1,2)) / (l(2,1)-l(1,1));
+ b = k(1,2) - a*k(1,1);
+ d = l(1,2) - c*l(1,1);
+
+ x = (d-b)/(a-c);
+ %y = a * x + b;
+ intersection = [x, a * x + b];
+end
+end
+
diff --git a/scripts/intersection_between_inclusive.m b/scripts/intersection_between_inclusive.m
new file mode 100644
index 0000000..eb54930
--- /dev/null
+++ b/scripts/intersection_between_inclusive.m
@@ -0,0 +1,59 @@
+function [ intersection ] = intersection_between_inclusive( k, vk, l, vl, between)
+%UNTITLED2 Calculates intersections of two lines specified by two sets of
+%points. If between is true, an empty value is returned if the intersection
+%is outside the endpoints of either one of the lines
+% lines in format [x1, y1; x2, y2]
+% [11, 12; 21, 22]
+% line k = ax + b, l = cx + d
+% vk = true iff k vertical
+
+if vk
+
+ c = (l(2,2) - l(1,2)) / (l(2,1)-l(1,1));
+ d = l(1,2) - c*l(1,1);
+ y = c * k(1) + d;
+ if isnan(y)
+ if vl
+ intersection = [];
+ else
+ error('non-parallel lines did not intersect');
+ end
+ elseif between && (y > max(k(3:4)) || y < min(k(3:4)) || ...
+ y > max(l(3:4)) || y < min(l(3:4)))
+ intersection = [];
+ else
+ intersection = [k(1), y];
+ end
+
+elseif vl
+ a = (k(2,2) - k(1,2)) / (k(2,1)-k(1,1));
+ b = k(1,2) - a*k(1,1);
+ x = l(1);
+ y = a * x + b;
+ if isnan(y)
+ error('non-parallel lines did not intersect');
+ elseif between && (y > max(l(3:4)) || y < (min(l(3:4))) || ...
+ y > max(k(3:4)) || y < (min(k(3:4))))
+ intersection = [];
+ else
+ intersection = [x, y];
+ end
+else
+ a = (k(2,2) - k(1,2)) / (k(2,1)-k(1,1));
+ c = (l(2,2) - l(1,2)) / (l(2,1)-l(1,1));
+ b = k(1,2) - a*k(1,1);
+ d = l(1,2) - c*l(1,1);
+
+ x = (d-b)/(a-c);
+ %y = a * x + b;
+ if isnan(x)
+ intersection = [];
+ elseif between && (x > max(l(1:2)) || x < (min(l(1:2))) || ...
+ x > max(k(1:2)) || x < (min(k(1:2))))
+ intersection = [];
+ else
+ intersection = [x, a * x + b];
+ end
+end
+end
+
diff --git a/scripts/intersections_of_cells.m b/scripts/intersections_of_cells.m
new file mode 100644
index 0000000..bcbce80
--- /dev/null
+++ b/scripts/intersections_of_cells.m
@@ -0,0 +1,40 @@
+function [ int_points ] = intersections_of_cells( line_points, ...
+ vertical_line, c_p_on, c_p_nw_on )
+%UNTITLED5 The matrix int_points contains the intersections with the line
+%specified by line_points for each of the cells defined by the points in
+%c_p_on. The format for the intersection points is x1,y1,x2,y2
+int_points = zeros(length(c_p_on),4);
+%length_on = zeros(length(c_p_on));
+for c = 1:length(c_p_on)
+ ind = 1; % places the points in the right position in the matrix line (c,:)
+ l = length(c_p_on{c});
+ for v = 1:l
+ if v < l
+ i = v+1; % could perhaps be handled with modulus, but =0 has to
+ %be avoided
+ else
+ i = 1;
+ end
+
+ % If the two points of a face are on opposite sides of the line,
+ % the face intersects the line:
+ if c_p_nw_on{c}(v) ~= c_p_nw_on{c}(i)
+ face_points = [c_p_on{c}(v,:); c_p_on{c}(i,:)];
+ vertical_face = abs(c_p_on{c}(v,1)- c_p_on{c}(i,1)) < 5*eps;
+ %find intersection points:
+ int_points(c,ind:ind+1) = intersection(line_points, ...
+ vertical_line, face_points, ...
+ vertical_face);
+
+ ind = ind + 2;
+ end
+ end
+
+end
+
+
+
+
+
+end
+
--
GitLab