changeset 939:46f5dc61d90b feature/utux2D

Update Schrodinger2d to work with new interface types, similar to LaplCurv.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 04 Dec 2018 14:25:03 -0800
parents 97291e1bd57c
children 31186559236d
files +scheme/Schrodinger2d.m
diffstat 1 files changed, 34 insertions(+), 25 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Schrodinger2d.m	Tue Dec 04 14:07:38 2018 -0800
+++ b/+scheme/Schrodinger2d.m	Tue Dec 04 14:25:03 2018 -0800
@@ -198,23 +198,25 @@
             end
         end
 
-        % opts     Struct that specifies the interface coupling.
+        % type     Struct that specifies the interface coupling.
         %          Fields:
-        %          -- interpolation:    struct of interpolation operators (empty for conforming grids)
-        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
-            if isempty(opts)
+        %          -- interpolation:    type of interpolation, default 'none'
+        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+
+            defaultType.interpolation = 'none';
+            default_struct('type', defaultType);
+
+            switch type.interpolation
+            case {'none', ''}
                 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary);
-            else
-                assertType(opts, 'struct');
-                if isfield(opts, 'I_local2neighbor')
-                    [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts);
-                else
-                    [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts);
-                end
+            case {'op','OP'}
+                [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
+            otherwise
+                error('Unknown type of interpolation: %s ', type.interpolation);
             end
         end
 
-        function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
+        function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
 
@@ -237,31 +239,38 @@
 
         end
 
-        function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
+        function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+
+            % User can request special interpolation operators by specifying type.interpOpSet
+            default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
+            interpOpSet = type.interpOpSet;
+
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
-
-            % Get boundary operators
-            [e_neighbour, d_neighbour] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
-            [e, d, H_gamma] = obj.get_boundary_ops(boundary);
+            [e_v, d_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
+            [e_u, d_u, H_gamma] = obj.get_boundary_ops(boundary);
             Hi = obj.Hi;
             a = obj.a;
 
             % Get outward unit normal component
             [~, n] = obj.get_boundary_number(boundary);
 
-            % Get interpolation operators from opts
-            I_u2v_good = opts.I_local2neighbor.good;
-            I_u2v_bad = opts.I_local2neighbor.bad;
-            I_v2u_good = opts.I_neighbor2local.good;
-            I_v2u_bad = opts.I_neighbor2local.bad;
+
+            % Find the number of grid points along the interface
+            m_u = size(e_u, 2);
+            m_v = size(e_v, 2);
+
+            % Build interpolation operators
+            intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
+            Iu2v = intOps.Iu2v;
+            Iv2u = intOps.Iv2u;
 
             sigma = -n*1i*a/2;
             tau = -n*(1i*a)'/2;
 
-            closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d';
-            penalty = -tau*Hi*d*H_gamma*I_v2u_good*e_neighbour' ...
-                      -sigma*Hi*e*H_gamma*I_v2u_bad*d_neighbour';
+            closure = tau*Hi*d_u*H_gamma*e_u' + sigma*Hi*e_u*H_gamma*d_u';
+            penalty = -tau*Hi*d_u*H_gamma*Iv2u.good*e_v' ...
+                      -sigma*Hi*e_u*H_gamma*Iv2u.bad*d_v';
 
         end