Mercurial > repos > public > sbplib

comparison +multiblock/DiffOp.m @ 498:324c927d8b1d feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
chnaged sbp interfacein 1d among many things
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 03 Mar 2017 16:19:41 +0100
parents 30ff8879162e
children 111fcbcff2e9
comparison
equal deleted inserted replaced
497:4905446f165e 498:324c927d8b1d
3 grid 3 grid
4 order 4 order
5 diffOps 5 diffOps
6 D 6 D
7 H 7 H
8 8
9 blockmatrixDiv 9 blockmatrixDiv
10 end 10 end
11 11
12 methods 12 methods
13 function obj = DiffOp(doHand, grid, order, doParam) 13 function obj = DiffOp(doHand, grid, order, doParam,timeDep)
14 % doHand -- may either be a function handle or a cell array of 14 % doHand -- may either be a function handle or a cell array of
15 % function handles for each grid. The function handle(s) 15 % function handles for each grid. The function handle(s)
16 % should be on the form do = doHand(grid, order, ...) 16 % should be on the form do = doHand(grid, order, ...)
17 % Additional parameters for each doHand may be provided in 17 % Additional parameters for each doHand may be provided in
18 % the doParam input. 18 % the doParam input.
23 % to the number of blocks then each element is sent to the 23 % to the number of blocks then each element is sent to the
24 % corresponding function handle as extra parameters: 24 % corresponding function handle as extra parameters:
25 % doHand(..., doParam{i}{:}) Otherwise doParam is sent as 25 % doHand(..., doParam{i}{:}) Otherwise doParam is sent as
26 % extra parameters to all doHand: doHand(..., doParam{:}) 26 % extra parameters to all doHand: doHand(..., doParam{:})
27 default_arg('doParam', []) 27 default_arg('doParam', [])
28 28
29 [getHand, getParam] = parseInput(doHand, grid, doParam); 29 [getHand, getParam] = parseInput(doHand, grid, doParam);
30 30
31 nBlocks = grid.nBlocks(); 31 nBlocks = grid.nBlocks();
32 32
33 obj.order = order; 33 obj.order = order;
34 34
35 % Create the diffOps for each block 35 % Create the diffOps for each block
36 obj.diffOps = cell(1, nBlocks); 36 obj.diffOps = cell(1, nBlocks);
37 for i = 1:nBlocks 37 for i = 1:nBlocks
38 h = getHand(i); 38 h = getHand(i);
39 p = getParam(i); 39 p = getParam(i);
40 if ~iscell(p) 40 if ~iscell(p)
41 p = {p}; 41 p = {p};
42 end 42 end
43 obj.diffOps{i} = h(grid.grids{i}, order, p{:}); 43 obj.diffOps{i} = h(grid.grids{i}, order, p{:});
44 end 44 end
45 45
46 46
47 % Build the norm matrix 47 % Build the norm matrix
48 H = cell(nBlocks, nBlocks); 48 H = cell(nBlocks, nBlocks);
49 for i = 1:nBlocks 49 for i = 1:nBlocks
50 H{i,i} = obj.diffOps{i}.H; 50 H{i,i} = obj.diffOps{i}.H;
51 end 51 end
52 obj.H = blockmatrix.toMatrix(H); 52 obj.H = blockmatrix.toMatrix(H);
53 53
54 54
55 % Build the differentiation matrix 55 % Build the differentiation matrix
56 obj.blockmatrixDiv = {grid.Ns, grid.Ns}; 56 switch timeDep
57 D = blockmatrix.zero(obj.blockmatrixDiv); 57 case {'n','no','N','No'}
58 for i = 1:nBlocks 58 obj.blockmatrixDiv = {grid.Ns, grid.Ns};
59 D{i,i} = obj.diffOps{i}.D; 59 D = blockmatrix.zero(obj.blockmatrixDiv);
60 end 60 for i = 1:nBlocks
61 61 D{i,i} = obj.diffOps{i}.D;
62 for i = 1:nBlocks 62 end
63 for j = 1:nBlocks 63
64 intf = grid.connections{i,j}; 64 for i = 1:nBlocks
65 if isempty(intf) 65 for j = 1:nBlocks
66 continue 66 intf = grid.connections{i,j};
67 end 67 if isempty(intf)
68 68 continue
69 69 end
70 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); 70
71 D{i,i} = D{i,i} + ii; 71
72 D{i,j} = D{i,j} + ij; 72 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2});
73 73 D{i,i} = D{i,i} + ii;
74 [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); 74 D{i,j} = D{i,j} + ij;
75 D{j,j} = D{j,j} + jj; 75
76 D{j,i} = D{j,i} + ji; 76 [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1});
77 end 77 D{j,j} = D{j,j} + jj;
78 end 78 D{j,i} = D{j,i} + ji;
79 obj.D = blockmatrix.toMatrix(D); 79 end
80 80 end
81 81 obj.D = blockmatrix.toMatrix(D);
82 case {'y','yes','Y','Yes'}
83 for i = 1:nBlocks
84 D{i,i} = @(t)obj.diffOps{i}.D(t);
85 end
86
87 for i = 1:nBlocks
88 for j = 1:nBlocks
89 intf = grid.connections{i,j};
90 if isempty(intf)
91 continue
92 end
93
94 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2});
95 D{i,i} = @(t)D{i,i}(t) + ii(t);
96 D{i,j} = ij;
97
98 [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1});
99 D{j,j} = @(t)D{j,j}(t) + jj(t);
100 D{j,i} = ji;
101 end
102 end
103 obj.D = D;
104 end
105
106
82 function [getHand, getParam] = parseInput(doHand, grid, doParam) 107 function [getHand, getParam] = parseInput(doHand, grid, doParam)
83 if ~isa(grid, 'multiblock.Grid') 108 if ~isa(grid, 'multiblock.Grid')
84 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); 109 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
85 end 110 end
86 111
87 if iscell(doHand) && length(doHand) == grid.nBlocks() 112 if iscell(doHand) && length(doHand) == grid.nBlocks()
88 getHand = @(i)doHand{i}; 113 getHand = @(i)doHand{i};
89 elseif isa(doHand, 'function_handle') 114 elseif isa(doHand, 'function_handle')
90 getHand = @(i)doHand; 115 getHand = @(i)doHand;
91 else 116 else
92 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); 117 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks');
93 end 118 end
94 119
95 if isempty(doParam) 120 if isempty(doParam)
96 getParam = @(i){}; 121 getParam = @(i){};
97 return 122 return
98 end 123 end
99 124
100 if ~iscell(doParam) 125 if ~iscell(doParam)
101 getParam = @(i)doParam; 126 getParam = @(i)doParam;
102 return 127 return
103 end 128 end
104 129
105 % doParam is a non-empty cell-array 130 % doParam is a non-empty cell-array
106 131
107 if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) 132 if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam))
108 % doParam is a cell-array of cell-arrays 133 % doParam is a cell-array of cell-arrays
109 getParam = @(i)doParam{i}; 134 getParam = @(i)doParam{i};
110 return 135 return
111 end 136 end
112 137
113 getParam = @(i)doParam; 138 getParam = @(i)doParam;
114 end 139 end
115 end 140 end
116 141
117 function ops = splitOp(obj, op) 142 function ops = splitOp(obj, op)
118 % Splits a matrix operator into a cell-matrix of matrix operators for 143 % Splits a matrix operator into a cell-matrix of matrix operators for
119 % each grid. 144 % each grid.
120 ops = sparse2cell(op, obj.NNN); 145 ops = sparse2cell(op, obj.NNN);
121 end 146 end
122 147
123 function op = getBoundaryOperator(obj, op, boundary) 148 function op = getBoundaryOperator(obj, op, boundary)
124 if iscell(boundary) 149 if iscell(boundary)
125 localOpName = [op '_' boundary{2}]; 150 localOpName = [op '_' boundary{2}];
126 blockId = boundary{1}; 151 blockId = boundary{1};
127 localOp = obj.diffOps{blockId}.(localOpName); 152 localOp = obj.diffOps{blockId}.(localOpName);
128 153
129 div = {obj.blockmatrixDiv{1}, size(localOp,2)}; 154 div = {obj.blockmatrixDiv{1}, size(localOp,2)};
130 blockOp = blockmatrix.zero(div); 155 blockOp = blockmatrix.zero(div);
131 blockOp{blockId,1} = localOp; 156 blockOp{blockId,1} = localOp;
132 op = blockmatrix.toMatrix(blockOp); 157 op = blockmatrix.toMatrix(blockOp);
133 return 158 return
134 else 159 else
135 % Boundary är en sträng med en boundary group i. 160 % Boundary är en sträng med en boundary group i.
136 end 161 end
137 end 162 end
138 163
139 % Creates the closure and penalty matrix for a given boundary condition, 164 % Creates the closure and penalty matrix for a given boundary condition,
140 % boundary -- the name of the boundary on the form {id,name} where 165 % boundary -- the name of the boundary on the form {id,name} where
141 % id is the number of a block and name is the name of a 166 % id is the number of a block and name is the name of a
142 % boundary of that block example: {1,'s'} or {3,'w'} 167 % boundary of that block example: {1,'s'} or {3,'w'}
143 function [closure, penalty] = boundary_condition(obj, boundary, type) 168 function [closure, penalty] = boundary_condition(obj, boundary, type)
144 I = boundary{1}; 169 I = boundary{1};
145 name = boundary{2}; 170 name = boundary{2};
146 171
147 % Get the closure and penaly matrices 172 % Get the closure and penaly matrices
148 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); 173 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type);
149 174
150 % Expand to matrix for full domain. 175 % Expand to matrix for full domain.
151 div = obj.blockmatrixDiv; 176 div = obj.blockmatrixDiv;
152 if ~iscell(blockClosure) 177 if ~iscell(blockClosure)
153 temp = blockmatrix.zero(div); 178 temp = blockmatrix.zero(div);
154 temp{I,I} = blockClosure; 179 temp{I,I} = blockClosure;
158 temp = blockmatrix.zero(div); 183 temp = blockmatrix.zero(div);
159 temp{I,I} = blockClosure{i}; 184 temp{I,I} = blockClosure{i};
160 closure{i} = blockmatrix.toMatrix(temp); 185 closure{i} = blockmatrix.toMatrix(temp);
161 end 186 end
162 end 187 end
163 188
164 div{2} = size(blockPenalty, 2); % Penalty is a column vector 189 div{2} = size(blockPenalty, 2); % Penalty is a column vector
165 if ~iscell(blockPenalty) 190 if ~iscell(blockPenalty)
166 p = blockmatrix.zero(div); 191 p = blockmatrix.zero(div);
167 p{I} = blockPenalty; 192 p{I} = blockPenalty;
168 penalty = blockmatrix.toMatrix(p); 193 penalty = blockmatrix.toMatrix(p);
172 p{I} = blockPenalty{i}; 197 p{I} = blockPenalty{i};
173 penalty{i} = blockmatrix.toMatrix(p); 198 penalty{i} = blockmatrix.toMatrix(p);
174 end 199 end
175 end 200 end
176 end 201 end
177 202
178 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 203 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
179 204
180 end 205 end
181 206
182 % Size returns the number of degrees of freedom 207 % Size returns the number of degrees of freedom
183 function N = size(obj) 208 function N = size(obj)
184 N = 0; 209 N = 0;
185 for i = 1:length(obj.diffOps) 210 for i = 1:length(obj.diffOps)
186 N = N + obj.diffOps{i}.size(); 211 N = N + obj.diffOps{i}.size();