Mercurial > repos > public > sbplib

comparison +multiblock/stitchSchemes.m @ 496:437fad4a47e1 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add 1d interface for shrodinger in 1d
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 24 Feb 2017 13:58:18 +0100
parents 419ec303e97d
children 324c927d8b1d
comparison
equal deleted inserted replaced
495:b91d23271481 496:437fad4a47e1
9 % each field with a parameter array that is sent to schm.boundary_condition 9 % each field with a parameter array that is sent to schm.boundary_condition
10 % 10 %
11 % Output parameters are cell arrays and cell matrices. 11 % Output parameters are cell arrays and cell matrices.
12 % 12 %
13 % Ex: [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound) 13 % Ex: [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound)
14 function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound) 14 function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound,timeDep)
15 default_arg('schmParam',[]); 15 default_arg('schmParam',[]);
16 16 default_arg('timeDep','N');
17 n_blocks = numel(grids); 17 n_blocks = numel(grids);
18 18
19 % Creating Schemes 19 % Creating Schemes
20 for i = 1:n_blocks 20 for i = 1:n_blocks
21 if isempty(schmParam); 21 if isempty(schmParam);
45 % Total norm 45 % Total norm
46 H = cell(n_blocks,n_blocks); 46 H = cell(n_blocks,n_blocks);
47 for i = 1:n_blocks 47 for i = 1:n_blocks
48 H{i,i} = schms{i}.H; 48 H{i,i} = schms{i}.H;
49 end 49 end
50 50
51 %% Total system matrix 51 %% Total system matrix
52 52
53 % Differentiation terms 53 % Differentiation terms
54 D = cell(n_blocks,n_blocks); 54 D = cell(n_blocks,n_blocks);
55 for i = 1:n_blocks 55 for i = 1:n_blocks
56 D{i,i} = schms{i}.D; 56 D{i,i} = schms{i}.D;
57 end 57 end
58 58
59 % Boundary penalty terms 59 % Boundary penalty terms
60 for i = 1:n_blocks 60 switch timeDep
61 if ~isstruct(bound{i}) 61 case {'n','no','N','No'}
62 continue 62 for i = 1:n_blocks
63 end 63 if ~isstruct(bound{i})
64 64 continue
65 fn = fieldnames(bound{i}); 65 end
66 for j = 1:length(fn); 66
67 bc = bound{i}.(fn{j}); 67 fn = fieldnames(bound{i});
68 if isempty(bc) 68 for j = 1:length(fn)
69 continue 69 bc = bound{i}.(fn{j});
70 if isempty(bc)
71 continue
72 end
73
74 [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:});
75 D{i,i} = D{i,i}+closure;
76 end
70 end 77 end
71 78
72 [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:}); 79 % Interface penalty terms
73 D{i,i} = D{i,i}+closure; 80 for i = 1:n_blocks
74 end 81 for j = 1:n_blocks
82 intf = conn{i,j};
83 if isempty(intf)
84 continue
85 end
86
87 [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2});
88 D{i,i} = D{i,i} + uu;
89 D{i,j} = uv;
90 D{j,j} = D{j,j} + vv;
91 D{j,i} = vu;
92 end
93 end
94 case {'y','yes','Y','Yes'}
95 for i = 1:n_blocks
96 if ~isstruct(bound{i})
97 continue
98 end
99
100 fn = fieldnames(bound{i});
101 for j = 1:length(fn)
102 bc = bound{i}.(fn{j});
103 if isempty(bc)
104 continue
105 end
106
107 [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:});
108 D{i,i} =@(t) D{i,i}(t) + closure(t);
109 end
110 end
111
112 % Interface penalty terms
113 for i = 1:n_blocks
114 for j = 1:n_blocks
115 intf = conn{i,j};
116 if isempty(intf)
117 continue
118 end
119
120 [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2});
121 D{i,i} = @(t) D{i,i}(t) + uu(t);
122 D{i,j} = uv;
123 D{j,j} = @(t)D{j,j}(t) + vv(t);
124 D{j,i} = vu;
125 end
126 end
75 end 127 end
76
77 % Interface penalty terms
78 for i = 1:n_blocks
79 for j = 1:n_blocks
80 intf = conn{i,j};
81 if isempty(intf)
82 continue
83 end
84
85 [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2});
86 D{i,i} = D{i,i} + uu;
87 D{i,j} = uv;
88 D{j,j} = D{j,j} + vv;
89 D{j,i} = vu;
90 end
91 end
92 end