Mercurial > repos > public > sbplib

comparison +multiblock/DiffOp.m @ 943:21394c78c72e feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Martin Almquist <malmquist@stanford.edu>
date Tue, 04 Dec 2018 15:24:36 -0800
parents 1c61d8fa9903 57760d7088ad
children a4ad90b37998 a3accd2f1283
comparison
equal deleted inserted replaced
942:35701c85e356 943:21394c78c72e
8 8
9 blockmatrixDiv 9 blockmatrixDiv
10 end 10 end
11 11
12 methods 12 methods
13 function obj = DiffOp(doHand, grid, order, doParam, intfTypes) 13 function obj = DiffOp(doHand, g, order, doParam, intfTypes)
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.
19 % grid -- a multiblock grid 19 % g -- a multiblock grid
20 % order -- integer specifying the order of accuracy 20 % order -- integer specifying the order of accuracy
21 % doParam -- may either be a cell array or a cell array of cell arrays 21 % doParam -- may either be a cell array or a cell array of cell arrays
22 % for each block. If it is a cell array with length equal 22 % for each block. If it is a cell array with length equal
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:
26 % extra parameters to all doHand: doHand(..., doParam{:}) 26 % extra parameters to all doHand: doHand(..., doParam{:})
27 % 27 %
28 % intfTypes (optional) -- nBlocks x nBlocks cell array of types for 28 % intfTypes (optional) -- nBlocks x nBlocks cell array of types for
29 % every interface. 29 % every interface.
30 default_arg('doParam', []) 30 default_arg('doParam', [])
31 default_arg('intfTypes', cell(grid.nBlocks(), grid.nBlocks()) ); 31 default_arg('intfTypes', cell(g.nBlocks(), g.nBlocks()) );
32 32
33 [getHand, getParam] = parseInput(doHand, grid, doParam); 33 [getHand, getParam] = parseInput(doHand, g, doParam);
34 34
35 obj.order = order; 35 obj.order = order;
36 nBlocks = grid.nBlocks(); 36 nBlocks = g.nBlocks();
37 37
38 % Create the diffOps for each block 38 % Create the diffOps for each block
39 obj.diffOps = cell(1, nBlocks); 39 obj.diffOps = cell(1, nBlocks);
40 for i = 1:nBlocks 40 for i = 1:nBlocks
41 h = getHand(i); 41 h = getHand(i);
42 p = getParam(i); 42 p = getParam(i);
43 if ~iscell(p) 43 if ~iscell(p)
44 p = {p}; 44 p = {p};
45 end 45 end
46 obj.diffOps{i} = h(grid.grids{i}, order, p{:}); 46 obj.diffOps{i} = h(g.grids{i}, order, p{:});
47 end 47 end
48 48
49 49
50 % Build the norm matrix 50 % Build the norm matrix
51 H = cell(nBlocks, nBlocks); 51 H = cell(nBlocks, nBlocks);
54 end 54 end
55 obj.H = blockmatrix.toMatrix(H); 55 obj.H = blockmatrix.toMatrix(H);
56 56
57 57
58 % Build the differentiation matrix 58 % Build the differentiation matrix
59 obj.blockmatrixDiv = {grid.Ns, grid.Ns}; 59 Ns = zeros(nBlocks,1);
60 for i = 1:nBlocks
61 Ns(i) = length(obj.diffOps{i}.D);
62 end
63 obj.blockmatrixDiv = {Ns, Ns};
60 D = blockmatrix.zero(obj.blockmatrixDiv); 64 D = blockmatrix.zero(obj.blockmatrixDiv);
61 for i = 1:nBlocks 65 for i = 1:nBlocks
62 D{i,i} = obj.diffOps{i}.D; 66 D{i,i} = obj.diffOps{i}.D;
63 end 67 end
64 68
65 for i = 1:nBlocks 69 for i = 1:nBlocks
66 for j = 1:nBlocks 70 for j = 1:nBlocks
67 intf = grid.connections{i,j}; 71 intf = g.connections{i,j};
68 if isempty(intf) 72 if isempty(intf)
69 continue 73 continue
70 end 74 end
71 75
72 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}, intfTypes{i,j}); 76 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}, intfTypes{i,j});
77 D{j,j} = D{j,j} + jj; 81 D{j,j} = D{j,j} + jj;
78 D{j,i} = D{j,i} + ji; 82 D{j,i} = D{j,i} + ji;
79 end 83 end
80 end 84 end
81 obj.D = blockmatrix.toMatrix(D); 85 obj.D = blockmatrix.toMatrix(D);
82 86 obj.grid = g;
83 87
84 function [getHand, getParam] = parseInput(doHand, grid, doParam) 88
85 if ~isa(grid, 'multiblock.Grid') 89 function [getHand, getParam] = parseInput(doHand, g, doParam)
90 if ~isa(g, 'multiblock.Grid')
86 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); 91 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
87 end 92 end
88 93
89 if iscell(doHand) && length(doHand) == grid.nBlocks() 94 if iscell(doHand) && length(doHand) == g.nBlocks()
90 getHand = @(i)doHand{i}; 95 getHand = @(i)doHand{i};
91 elseif isa(doHand, 'function_handle') 96 elseif isa(doHand, 'function_handle')
92 getHand = @(i)doHand; 97 getHand = @(i)doHand;
93 else 98 else
94 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); 99 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks');
104 return 109 return
105 end 110 end
106 111
107 % doParam is a non-empty cell-array 112 % doParam is a non-empty cell-array
108 113
109 if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) 114 if length(doParam) == g.nBlocks() && all(cellfun(@iscell, doParam))
110 % doParam is a cell-array of cell-arrays 115 % doParam is a cell-array of cell-arrays
111 getParam = @(i)doParam{i}; 116 getParam = @(i)doParam{i};
112 return 117 return
113 end 118 end
114 119
138 case 'multiblock.BoundaryGroup' 143 case 'multiblock.BoundaryGroup'
139 op = sparse(size(obj.D,1),0); 144 op = sparse(size(obj.D,1),0);
140 for i = 1:length(boundary) 145 for i = 1:length(boundary)
141 op = [op, obj.getBoundaryOperator(opName, boundary{i})]; 146 op = [op, obj.getBoundaryOperator(opName, boundary{i})];
142 end 147 end
148 otherwise
149 error('Unknown boundary indentifier')
150 end
151 end
152
153 function op = getBoundaryQuadrature(obj, boundary)
154 opName = 'H';
155 switch class(boundary)
156 case 'cell'
157 localOpName = [opName '_' boundary{2}];
158 blockId = boundary{1};
159 op = obj.diffOps{blockId}.(localOpName);
160
161 return
162 case 'multiblock.BoundaryGroup'
163 N = length(boundary);
164 H_bm = cell(N,N);
165 for i = 1:N
166 H_bm{i,i} = obj.getBoundaryQuadrature(boundary{i});
167 end
168 op = blockmatrix.toMatrix(H_bm);
143 otherwise 169 otherwise
144 error('Unknown boundary indentifier') 170 error('Unknown boundary indentifier')
145 end 171 end
146 end 172 end
147 173
175 201
176 % Get the closure and penaly matrices 202 % Get the closure and penaly matrices
177 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); 203 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type);
178 204
179 % Expand to matrix for full domain. 205 % Expand to matrix for full domain.
180 div = obj.blockmatrixDiv; 206 closure = multiblock.local2globalClosure(blockClosure, obj.blockmatrixDiv, I);
181 if ~iscell(blockClosure) 207 penalty = multiblock.local2globalPenalty(blockPenalty, obj.blockmatrixDiv, I);
182 temp = blockmatrix.zero(div);
183 temp{I,I} = blockClosure;
184 closure = blockmatrix.toMatrix(temp);
185 else
186 for i = 1:length(blockClosure)
187 temp = blockmatrix.zero(div);
188 temp{I,I} = blockClosure{i};
189 closure{i} = blockmatrix.toMatrix(temp);
190 end
191 end
192
193 if ~iscell(blockPenalty)
194 div{2} = size(blockPenalty, 2); % Penalty is a column vector
195 p = blockmatrix.zero(div);
196 p{I} = blockPenalty;
197 penalty = blockmatrix.toMatrix(p);
198 else
199 % TODO: used by beam equation, should be eliminated. SHould only set one BC per call
200 for i = 1:length(blockPenalty)
201 div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector
202 p = blockmatrix.zero(div);
203 p{I} = blockPenalty{i};
204 penalty{i} = blockmatrix.toMatrix(p);
205 end
206 end
207 end 208 end
208 209
209 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 210 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
210 error('not implemented') 211 error('not implemented')
211 end 212 end