Upload new file

scripts/line_norm.m

+function [relativeError, errorSquared, normalizationSquared] = ...
+                                    line_norm(fracXAndP, p, t, x, FRAC)
+%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.
+
+boundary_line = false;
+
+
+if fracXAndP(1)==fracXAndP(2)
+    endPoints  = [fracXAndP(1),min(fracXAndP(:,2)); ...
+                  fracXAndP(1),max(fracXAndP(:,2))];
+    a=endPoints(1);
+    b=0;
+    isVertical = true;
+    
+    fractureLength = endPoints(4)-endPoints(3);
+else
+    [x1,iMin]       = min(fracXAndP(:,1));
+    [x2,iMax]       = max(fracXAndP(:,1));
+    endPoints       = [x1,fracXAndP(iMin,2);x2,fracXAndP(iMax,2)];
+    a               = (endPoints(2,2) - endPoints(1,2)) / (endPoints(2,1)-endPoints(1,1));
+    b               = endPoints(1,2) - a*endPoints(1,1);
+    isVertical      = false;
+    fractureLength  = hypot(abs(endPoints(1)-endPoints(2)), abs(endPoints(3)-endPoints(4)));
+end
+
+[verticesNorthWest, verticesOn] = check_points(x,a,b,isVertical);
+
+if ~FRAC
+    [VerticesNorthWest, VerticesOn] = check_points(X,a,b,isVertical);
+end   
+
+
+% The code must unfortunately allow for cells with different number of vertices 
+% both for fine grid and coarse for generality. 
+% Therefore, cellPointsOn is a CELL with nVertices(c) pointer extracted
+% from t. cellPointsNorthWestOn is the corresponding logical northWest
+% cell, true for points north west of the fracture.
+% Triangulation faces might coincide with the fracture:
+[cellPointsOn, pOn, cellPointsNorthWestOn, cellVerticesOn, verticesP] = ... 
+        find_on_line(t, x, verticesNorthWest, verticesOn, p);
+if ~FRAC
+    [CellPointsOn, POn, CellPointsNorthWestOn, CellVerticesOn, VerticesP] = ...
+            find_on_line(T, X, VerticesNorthWest, VerticesOn, P);
+else
+    CellVerticesOn = [];
+    VerticesP = [];
+end
+
+% Those faces will not be unique. The values corresponding to the two 
+% neighbouring cells are merged:
+if ~boundary_line && ~isempty(cellVerticesOn)
+    [cellVerticesOn, verticesP] = combine_faces(cellVerticesOn, verticesP,isVertical);
+end
+if ~boundary_line && ~isempty(CellVerticesOn)    
+    [CellVerticesOn, VerticesP] = combine_faces(CellVerticesOn, VerticesP,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, ...
+                                    cellPointsOn, cellPointsNorthWestOn);
+                                
+if FRAC
+    nCells = size(fracXAndP,1);
+    ind_two = linspace(2,nCells,nCells/2);
+    ind_one = ind_two-1;
+    IntersectionPoints = [fracXAndP(ind_one,1:2), fracXAndP(ind_two,1:2)];
+    POn = fracXAndP(ind_two,3);
+else
+    [IntersectionPoints] = intersections_of_cells(endPoints, isVertical, ...
+                                    CellPointsOn, CellPointsNorthWestOn);
+end        
+
+% Add the values for the faces coinciding with the line:
+if ~isempty(cellVerticesOn)
+    intersectionPoints = [intersectionPoints;cellVerticesOn];
+    pOn = [pOn;verticesP];
+end
+if ~isempty(CellVerticesOn)
+    IntersectionPoints = [IntersectionPoints;CellVerticesOn];
+    POn = [POn;VerticesP];
+end
+
+
+% Loop through large cells, find smalls cells (partly) inside and evaluate
+% norms in evaluate_norm:
+errorSquared = evaluate_norm(IntersectionPoints, intersectionPoints, POn, pOn, isVertical,endPoints);
+
+normalizationSquared = (max(p)-min(p))^2*fractureLength;
+relativeError = sqrt((errorSquared)/(normalizationSquared));
+
+
+end
+