Mercurial > repos > public > sbplib

changeset 998:2b1b944deae1 feature/getBoundaryOp

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add getBoundaryOperator to all 1d schemes. Did not add getBoundaryQuadrature because it doesnt make sense in 1d (?)
author Martin Almquist <malmquist@stanford.edu>
date Sat, 12 Jan 2019 13:35:19 -0800
parents 78db023a7fe3
children 337c4d1dcef5
files +scheme/Beam.m +scheme/Euler1d.m +scheme/Laplace1d.m +scheme/Utux.m
diffstat 4 files changed, 195 insertions(+), 50 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Beam.m	Sat Jan 12 11:57:50 2019 -0800
+++ b/+scheme/Beam.m	Sat Jan 12 13:35:19 2019 -0800
@@ -86,7 +86,8 @@
         function [closure, penalty] = boundary_condition(obj,boundary,type)
             default_arg('type','dn');
 
-            [e, d1, d2, d3, s] = obj.get_boundary_ops(boundary);
+            [e, d1, d2, d3] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary);
+            s = obj.getBoundarySign(boundary);
             gamm = obj.gamm;
             delt = obj.delt;
 
@@ -173,14 +174,15 @@
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type)
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
-            [e_u,d1_u,d2_u,d3_u,s_u] = obj.get_boundary_ops(boundary);
-            [e_v,d1_v,d2_v,d3_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
+            [e_u, d1_u, d2_u, d3_u] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary);
+            s_u = obj.getBoundarySign(boundary);
 
+            [e_v, d1_v, d2_v, d3_v] = neighbour_scheme.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, neighbour_boundary);
+            s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
 
             alpha_u = obj.alpha;
             alpha_v = neighbour_scheme.alpha;
 
-
             switch boundary
                 case 'l'
                     interface_opt = obj.opt.interface_l;
@@ -234,22 +236,70 @@
             penalty = -obj.Hi*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v');
         end
 
-        % Returns the boundary ops and sign for the boundary specified by the string boundary.
-        % The right boundary is considered the positive boundary
-        function [e, d1, d2, d3, s] = get_boundary_ops(obj,boundary)
+        % Returns the boundary operator op for the boundary specified by the string boundary.
+        % op        -- string or a cell array of strings
+        % boundary  -- string
+        function varargout = getBoundaryOperator(obj, op, boundary)
+
+            if ~ismember(boundary, {'l', 'r'})
+                error('No such boundary: boundary = %s',boundary);
+            end
+
+            if ~iscell(op)
+                op = {op};
+            end
+
+            for i = 1:numel(op)
+                switch op{i}
+                case 'e'
+                    switch boundary
+                    case 'l'
+                        e = obj.e_l;
+                    case 'r'
+                        e = obj.e_r;
+                    end
+                    varargout{i} = e;
+
+                case 'd1'
+                    switch boundary
+                    case 'l'
+                        d1 = obj.d1_l;
+                    case 'r'
+                        d1 = obj.d1_r;
+                    end
+                    varargout{i} = d1;
+                end
+
+                case 'd2'
+                    switch boundary
+                    case 'l'
+                        d2 = obj.d2_l;
+                    case 'r'
+                        d2 = obj.d2_r;
+                    end
+                    varargout{i} = d2;
+                end
+
+                case 'd3'
+                    switch boundary
+                    case 'l'
+                        d3 = obj.d3_l;
+                    case 'r'
+                        d3 = obj.d3_r;
+                    end
+                    varargout{i} = d3;
+                end
+            end
+        end
+
+        % Returns the boundary sign. The right boundary is considered the positive boundary
+        % boundary -- string
+        function s = getBoundarySign(obj, boundary)
             switch boundary
-                case 'l'
-                    e  = obj.e_l;
-                    d1 = obj.d1_l;
-                    d2 = obj.d2_l;
-                    d3 = obj.d3_l;
+                case {'r'}
+                    s = 1;
+                case {'l'}
                     s = -1;
-                case 'r'
-                    e  = obj.e_r;
-                    d1 = obj.d1_r;
-                    d2 = obj.d2_r;
-                    d3 = obj.d3_r;
-                    s = 1;
                 otherwise
                     error('No such boundary: boundary = %s',boundary);
             end
--- a/+scheme/Euler1d.m	Sat Jan 12 11:57:50 2019 -0800
+++ b/+scheme/Euler1d.m	Sat Jan 12 13:35:19 2019 -0800
@@ -201,7 +201,8 @@
         % Enforces the boundary conditions
         %  w+ = R*w- + g(t)
         function closure = boundary_condition(obj,boundary, type, varargin)
-            [e_s,e_S,s] = obj.get_boundary_ops(boundary);
+            [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary);
+            s = obj.getBoundarySign(boundary);
 
             % Boundary condition on form
             %   w_in = R*w_out + g,       where g is data
@@ -232,7 +233,8 @@
         %
         % Returns closure(q,g)
         function closure = boundary_condition_L(obj, boundary, L_fun, p_in)
-            [e_s,e_S,s] = obj.get_boundary_ops(boundary);
+            [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary);
+            s = obj.getBoundarySign(boundary);
 
             p_ot = 1:3;
             p_ot(p_in) = [];
@@ -273,7 +275,8 @@
 
         % Return closure(q,g)
         function closure = boundary_condition_char(obj,boundary)
-            [e_s,e_S,s] = obj.get_boundary_ops(boundary);
+            [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary);
+            s = obj.getBoundarySign(boundary);
 
             function o = closure_fun(q, w_data)
                 q_s = e_S'*q;
@@ -314,7 +317,7 @@
 
         % Return closure(q,[v; p])
         function closure = boundary_condition_inflow(obj, boundary)
-            [~,~,s] = obj.get_boundary_ops(boundary);
+            s = obj.getBoundarySign(boundary);
 
              switch s
                 case -1
@@ -335,7 +338,7 @@
 
         % Return closure(q, p)
         function closure = boundary_condition_outflow(obj, boundary)
-            [~,~,s] = obj.get_boundary_ops(boundary);
+            s = obj.getBoundarySign(boundary);
 
              switch s
                 case -1
@@ -352,7 +355,7 @@
 
         % Return closure(q,[v; rho])
         function closure = boundary_condition_inflow_rho(obj, boundary)
-            [~,~,s] = obj.get_boundary_ops(boundary);
+            s = obj.getBoundarySign(boundary);
 
              switch s
                 case -1
@@ -372,7 +375,7 @@
 
         % Return closure(q,rho)
         function closure = boundary_condition_outflow_rho(obj, boundary)
-            [~,~,s] = obj.get_boundary_ops(boundary);
+            s = obj.getBoundarySign(boundary);
 
              switch s
                 case -1
@@ -388,7 +391,8 @@
 
         % Set wall boundary condition v = 0.
         function closure = boundary_condition_wall(obj,boundary)
-            [e_s,e_S,s] = obj.get_boundary_ops(boundary);
+            [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary);
+            s = obj.getBoundarySign(boundary);
 
             % Vill vi sätta penalty på karateristikan som är nära noll också eller vill
             % vi låta den vara fri?
@@ -478,18 +482,50 @@
             penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v');
         end
 
-        % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
-        % The right boundary is considered the positive boundary
-        function [e,E,s] = get_boundary_ops(obj,boundary)
+        % Returns the boundary operator op for the boundary specified by the string boundary.
+        % op        -- string or a cell array of strings
+        % boundary  -- string
+        function varargout = getBoundaryOperator(obj, op, boundary)
+
+            if ~iscell(op)
+                op = {op};
+            end
+
+            for i = 1:numel(op)
+                switch op{i}
+                case 'e'
+                    switch boundary
+                    case 'l'
+                        e = obj.e_l;
+                    case 'r'
+                        e = obj.e_r;
+                    otherwise
+                        error('No such boundary: boundary = %s',boundary);
+                    end
+                    varargout{i} = e;
+
+                case 'E'
+                    switch boundary
+                    case 'l'
+                        E = obj.e_L;
+                    case 'r'
+                        E = obj.e_R;
+                    otherwise
+                        error('No such boundary: boundary = %s',boundary);
+                    end
+                    varargout{i} = E;
+                end
+            end
+        end
+
+        % Returns the boundary sign. The right boundary is considered the positive boundary
+        % boundary -- string
+        function s = getBoundarySign(obj, boundary)
             switch boundary
-                case 'l'
-                    e  = obj.e_l;
-                    E  = obj.e_L;
+                case {'r'}
+                    s = 1;
+                case {'l'}
                     s = -1;
-                case 'r'
-                    e  = obj.e_r;
-                    E  = obj.e_R;
-                    s = 1;
                 otherwise
                     error('No such boundary: boundary = %s',boundary);
             end
--- a/+scheme/Laplace1d.m	Sat Jan 12 11:57:50 2019 -0800
+++ b/+scheme/Laplace1d.m	Sat Jan 12 13:35:19 2019 -0800
@@ -56,7 +56,8 @@
             default_arg('type','neumann');
             default_arg('data',0);
 
-            [e,d,s] = obj.get_boundary_ops(boundary);
+            [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
+            s = obj.getBoundarySign(boundary);
 
             switch type
                 % Dirichlet boundary condition
@@ -86,10 +87,11 @@
         function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
+            [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
+            s_u = obj.getBoundarySign(boundary);
 
-            [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
-            [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
-
+            [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
+            s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
 
             a_u = obj.a;
             a_v = neighbour_scheme.a;
@@ -111,18 +113,50 @@
             penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v');
         end
 
-        % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
-        % The right boundary is considered the positive boundary
-        function [e,d,s] = get_boundary_ops(obj,boundary)
+        % Returns the boundary operator op for the boundary specified by the string boundary.
+        % op        -- string or a cell array of strings
+        % boundary  -- string
+        function varargout = getBoundaryOperator(obj, op, boundary)
+
+            if ~iscell(op)
+                op = {op};
+            end
+
+            for i = 1:numel(op)
+                switch op{i}
+                case 'e'
+                    switch boundary
+                    case 'l'
+                        e = obj.e_l;
+                    case 'r'
+                        e = obj.e_r;
+                    otherwise
+                        error('No such boundary: boundary = %s',boundary);
+                    end
+                    varargout{i} = e;
+
+                case 'd'
+                    switch boundary
+                    case 'l'
+                        d = obj.d_l;
+                    case 'r'
+                        d = obj.d_r;
+                    otherwise
+                        error('No such boundary: boundary = %s',boundary);
+                    end
+                    varargout{i} = d;
+                end
+            end
+        end
+
+        % Returns the boundary sign. The right boundary is considered the positive boundary
+        % boundary -- string
+        function s = getBoundarySign(obj, boundary)
             switch boundary
-                case 'l'
-                    e = obj.e_l;
-                    d = obj.d_l;
+                case {'r'}
+                    s = 1;
+                case {'l'}
                     s = -1;
-                case 'r'
-                    e = obj.e_r;
-                    d = obj.d_r;
-                    s = 1;
                 otherwise
                     error('No such boundary: boundary = %s',boundary);
             end
--- a/+scheme/Utux.m	Sat Jan 12 11:57:50 2019 -0800
+++ b/+scheme/Utux.m	Sat Jan 12 13:35:19 2019 -0800
@@ -72,6 +72,31 @@
 
          end
 
+        % Returns the boundary operator op for the boundary specified by the string boundary.
+        % op        -- string or a cell array of strings
+        % boundary  -- string
+        function varargout = getBoundaryOperator(obj, op, boundary)
+
+            if ~iscell(op)
+                op = {op};
+            end
+
+            for i = 1:numel(op)
+                switch op{i}
+                case 'e'
+                    switch boundary
+                    case 'l'
+                        e = obj.e_l;
+                    case 'r'
+                        e = obj.e_r;
+                    otherwise
+                        error('No such boundary: boundary = %s',boundary);
+                    end
+                    varargout{i} = e;
+                end
+            end
+        end
+
         function N = size(obj)
             N = obj.m;
         end