Replace evaluate_norm.m

scripts/evaluate_norm.m

-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