Mercurial > repos > public > sbplib
changeset 1139:6bc93c091682 feature/laplace_curvilinear_test
Implement minimum check in new scheme.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Fri, 21 Jun 2019 16:27:49 +0200 |
| parents | afd06a84b69c |
| children | 738de3a4058b |
| files | +scheme/LaplaceCurvilinearMin.m |
| diffstat | 1 files changed, 117 insertions(+), 14 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/LaplaceCurvilinearMin.m Mon Jun 10 14:39:14 2019 +0200 +++ b/+scheme/LaplaceCurvilinearMin.m Fri Jun 21 16:27:49 2019 +0200 @@ -47,6 +47,9 @@ % Temporary, only used for nonconforming interfaces but should be removed. lambda + + % Number of boundary points in minumum check + bp end methods @@ -58,6 +61,16 @@ default_arg('a', 1); default_arg('b', @(x,y) 0*x + 1); + % Number of boundary points in minimum check + switch order + case 2 + obj.bp = 2; + case 4 + obj.bp = 4; + case 6 + obj.bp = 7; + end + % assert(isa(g, 'grid.Curvilinear')) if isa(a, 'function_handle') a = grid.evalOn(g, a); @@ -285,24 +298,54 @@ Hi = obj.Hi; a = obj.a; b_b = e'*obj.b*e; + b = obj.b; switch type % Dirichlet boundary condition case {'D','d','dirichlet'} tuning = 1.0; - sigma = 0*b_b; + sigma = 0*b; for i = 1:obj.dim - sigma = sigma + e'*J*K{i,m}*K{i,m}*e; + sigma = sigma + b*J*K{i,m}*K{i,m}; end + + % Minimum check on sigma + mx = obj.m(1); + my = obj.m(2); + bp = obj.bp; + sigma_mat = reshape(diag(sigma), my, mx); + switch boundary + case 'w' + sigma_min = min(sigma_mat(:,1:bp), [], 2); + sigma_mat = repmat(sigma_min, 1, mx); + case 'e' + sigma_min = min(sigma_mat(:,end-bp+1:end), [], 2); + sigma_mat = repmat(sigma_min, 1, mx); + case 's' + sigma_min = min(sigma_mat(1:bp,:), [], 1); + sigma_mat = repmat(sigma_min, my, 1); + case 'n' + sigma_min = min(sigma_mat(end-bp+1:end,:), [], 1); + sigma_mat = repmat(sigma_min, my, 1); + end + sigma_min = sigma_mat(:); + sigma_min = spdiag(sigma_min); + + % Window + sigma = e'*sigma*e; + sigma_min = e'*sigma_min*e; + sigma = sigma/s_b; - tau = tuning*(1/th_R + obj.dim/th_H)*sigma; + sigma_min = sigma_min/s_b; + + tau = tuning*(1/th_R*sigma/sigma_min*sigma + obj.dim/th_H*sigma); closure = a*Hi*d*b_b*H_b*e' ... - -a*Hi*e*tau*b_b*H_b*e'; + -a*Hi*e*tau*H_b*e'; penalty = -a*Hi*d*b_b*H_b ... - +a*Hi*e*tau*b_b*H_b; + +a*Hi*e*tau*H_b; % Neumann boundary condition. Note that the penalty is for du/dn and not b*du/dn. @@ -362,6 +405,7 @@ % Coefficients, u K_u = u.K; J_u = u.J; + b_u = u.b; b_b_u = e_u'*u.b*e_u; % Boundary operators, v @@ -381,25 +425,84 @@ % Coefficients, v K_v = v.K; J_v = v.J; + b_v = v.b; b_b_v = e_v'*v.b*e_v; %--- Penalty strength tau ------------- - sigma_u = 0*b_b_u; - sigma_v = 0*b_b_v; + sigma_u = 0*b_u; + sigma_v = 0*b_v; for i = 1:obj.dim - sigma_u = sigma_u + e_u'*J_u*K_u{i,m_u}*K_u{i,m_u}*e_u; - sigma_v = sigma_v + e_v'*J_v*K_v{i,m_v}*K_v{i,m_v}*e_v; + sigma_u = sigma_u + b_u*J_u*K_u{i,m_u}*K_u{i,m_u}; + sigma_v = sigma_v + b_v*J_v*K_v{i,m_v}*K_v{i,m_v}; + end + + %--- Minimum check on sigma_u ---- + mx = u.m(1); + my = u.m(2); + bp = u.bp; + sigma_mat = reshape(diag(sigma_u), my, mx); + switch boundary + case 'w' + sigma_min = min(sigma_mat(:,1:bp), [], 2); + sigma_mat = repmat(sigma_min, 1, mx); + case 'e' + sigma_min = min(sigma_mat(:,end-bp+1:end), [], 2); + sigma_mat = repmat(sigma_min, 1, mx); + case 's' + sigma_min = min(sigma_mat(1:bp,:), [], 1); + sigma_mat = repmat(sigma_min, my, 1); + case 'n' + sigma_min = min(sigma_mat(end-bp+1:end,:), [], 1); + sigma_mat = repmat(sigma_min, my, 1); end - sigma_u = sigma_u/s_b_u; - sigma_v = sigma_v/s_b_v; + sigma_min_u = sigma_mat(:); + sigma_min_u = spdiag(sigma_min_u); + + % Window + sigma_u = e_u'*sigma_u*e_u; + sigma_min_u = e_u'*sigma_min_u*e_u; + % ------------------------------- - tau_R_u = 1/th_R_u*sigma_u; - tau_R_v = 1/th_R_v*sigma_v; + %--- Minimum check on sigma_v ---- + mx = v.m(1); + my = v.m(2); + bp = v.bp; + sigma_mat = reshape(diag(sigma_v), my, mx); + switch neighbour_boundary + case 'w' + sigma_min = min(sigma_mat(:,1:bp), [], 2); + sigma_mat = repmat(sigma_min, 1, mx); + case 'e' + sigma_min = min(sigma_mat(:,end-bp+1:end), [], 2); + sigma_mat = repmat(sigma_min, 1, mx); + case 's' + sigma_min = min(sigma_mat(1:bp,:), [], 1); + sigma_mat = repmat(sigma_min, my, 1); + case 'n' + sigma_min = min(sigma_mat(end-bp+1:end,:), [], 1); + sigma_mat = repmat(sigma_min, my, 1); + end + sigma_min_v = sigma_mat(:); + sigma_min_v = spdiag(sigma_min_v); + + % Window + sigma_v = e_v'*sigma_v*e_v; + sigma_min_v = e_v'*sigma_min_v*e_v; + % ------------------------------- + + sigma_u = sigma_u/s_b_u; + sigma_min_u = sigma_min_u/s_b_u; + + sigma_v = sigma_v/s_b_v; + sigma_min_v = sigma_min_v/s_b_v; + + tau_R_u = 1/th_R_u*sigma_u/sigma_min_u*sigma_u; + tau_R_v = 1/th_R_v*sigma_v/sigma_min_v*sigma_v; tau_H_u = dim*1/th_H_u*sigma_u; tau_H_v = dim*1/th_H_v*sigma_v; - tau = 1/4*tuning*(b_b_u*(tau_R_u + tau_H_u) + b_b_v*(tau_R_v + tau_H_v)); + tau = 1/4*tuning*((tau_R_u + tau_H_u) + (tau_R_v + tau_H_v)); %-------------------------------------- % Operators/coefficients that are only required from this side