From a72676a706abfc1eb710b828314957ceb6d58c5f Mon Sep 17 00:00:00 2001 From: Ivar Stefansson <ist050@uib.no> Date: Wed, 21 Dec 2016 18:29:34 +0100 Subject: [PATCH] Replace evaluate_norm.m --- scripts/evaluate_norm.m | 383 ++++++++++++++++++++-------------------- 1 file changed, 192 insertions(+), 191 deletions(-) diff --git a/scripts/evaluate_norm.m b/scripts/evaluate_norm.m index 15be47b..db486da 100644 --- a/scripts/evaluate_norm.m +++ b/scripts/evaluate_norm.m @@ -1,192 +1,193 @@ -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 +function [ errorSquared ] = evaluate_norm( IntersectionPoints, intersectionPoints, ... + P, p, isVertical,endPoints) +%evaluate_norm returns L2norm^2 of the pressure solutions along a line. +% The points IntersectionPoints and intersectionPoints all lie on the line. +% The length of the line is determined by endPoints. + +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 endPoints(1) ~= endPoints(2) + if isVertical + first = find(ascending(:,2) >= endPoints(3),1); + last = find(ascending(:,4) <= endPoints(4),1,'last'); + ascending = ascending(first:last,:); + + First = find(Ascending(:,2) >= endPoints(3),1); + Last = find(Ascending(:,4) <= endPoints(4),1,'last'); + Ascending = Ascending(First:Last,:); + else + first = find(ascending(:,1) >= endPoints(1),1); + last = find(ascending(:,3) <= endPoints(2),1,'last'); + ascending = ascending(first:last,:); + + First = find(Ascending(:,1) >= endPoints(1),1); + Last = find(Ascending(:,3) <= endPoints(2),1,'last'); + Ascending = Ascending(First:Last,:); + + end +end +intersectionPoints = ascending(:,1:4); +p = ascending(:,5); +IntersectionPoints = Ascending(:,1:4); +P = Ascending(:,5); + +errorSquared = 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); + areIn = 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(areIn) + + + % fine cell is larger (on the line) than coarse cell + elseif length(areIn) == 1 + largerLength = Yc(2)-Yc(1); + errorSquared = errorSquared + largerLength*(p(areIn)-P(C))^2; + + % at least two fine cells (partly) inside the coarse + else + bottom = areIn(1); + top = areIn(end); + bottomLength = intersectionPoints(bottom,4) - Yc(1); + topLength = Yc(2) - intersectionPoints(top,2); + errorSquared = errorSquared + sum([bottomLength,topLength].* ... + ([p(bottom),p(top)]-P(C)).^2); + + completelyIn = areIn(2:end-1); + if ~isempty(completelyIn) + completelyLength = intersectionPoints(completelyIn, 4) ... + -intersectionPoints(completelyIn,2); + errorSquared = errorSquared + sum(completelyLength.*(p(completelyIn)-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); + areIn = 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(areIn) + + + % fine cell is larger (on the line) than coarse cell + elseif length(areIn) == 1 + largerLength = Xc(2)-Xc(1); + errorSquared = errorSquared + largerLength*(p(areIn)-P(C))^2; + % at least two fine cells (partly) inside the coarse + else + bottom = areIn(1); + top = areIn(end); + bottomLength = intersectionPoints(bottom,3) - Xc(1); + topLength = Xc(2) - intersectionPoints(top,1) ; + errorSquared = errorSquared + sum([bottomLength,topLength].* ... + ([p(bottom),p(top)]-P(C)).^2); + + + completelyIn = areIn(2:end-1); + if ~isempty(completelyIn) + completelyLength = intersectionPoints(completelyIn, 3) ... + -intersectionPoints(completelyIn,1); + errorSquared = errorSquared + sum(completelyLength.*(p(completelyIn)-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); + areIn = 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(areIn) + + + % fine cell is larger (on the line) than coarse cell + elseif length(areIn) == 1 + largerLength = norm(IntersectionPoints(C,1:2) - ... + IntersectionPoints(C,3:4), 2); + errorSquared = errorSquared + largerLength*(p(areIn)-P(C))^2; + % at least two fine cells (partly) inside the coarse + else + bottom = areIn(1); + top = areIn(end); + + bottomLength = norm(intersectionPoints(bottom,3:4) - ... + [Xc(1), Yc(1)], 2); + topLength = norm(intersectionPoints(top,1:2) - ... + [Xc(2), Yc(2)], 2); + errorSquared = errorSquared + sum([bottomLength,topLength].* ... + ([p(bottom),p(top)]-P(C)).^2); + + completelyIn = areIn(2:end-1); + if ~isempty(completelyIn) + completelyLength = hypot(intersectionPoints(completelyIn, 1) ... + -intersectionPoints(completelyIn,3), ... + intersectionPoints(completelyIn, 2) - ... + intersectionPoints(completelyIn,4)); + errorSquared = errorSquared + sum(completelyLength.*(p(completelyIn)-P(C)).^2); + end + end + end end \ No newline at end of file -- GitLab