diff --git a/scripts/Exact2dNorm.m b/scripts/Exact2dNorm.m deleted file mode 100644 index 5c89cc5867acec0f4803e38d0d13d44e3c5cd1bd..0000000000000000000000000000000000000000 --- a/scripts/Exact2dNorm.m +++ /dev/null @@ -1,117 +0,0 @@ -%% 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')