Mercurial > repos > public > sbplib

comparison +multiblock/DiffOp.m @ 1072:6468a5f6ec79 feature/grids/LaplaceSquared

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 12 Feb 2019 17:12:42 +0100
parents 9c8ed00732fd
children 60c875c18de3
comparison
equal deleted inserted replaced
1071:92cb03e64ca4 1072:6468a5f6ec79
8 8
9 blockmatrixDiv 9 blockmatrixDiv
10 end 10 end
11 11
12 methods 12 methods
13 function obj = DiffOp(doHand, g, order, doParam) 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.
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:
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 %
28 % intfTypes (optional) -- nBlocks x nBlocks cell array of types for
29 % every interface.
27 default_arg('doParam', []) 30 default_arg('doParam', [])
31 default_arg('intfTypes', cell(g.nBlocks(), g.nBlocks()) );
28 32
29 [getHand, getParam] = parseInput(doHand, g, doParam); 33 [getHand, getParam] = parseInput(doHand, g, doParam);
30 34
35 obj.order = order;
31 nBlocks = g.nBlocks(); 36 nBlocks = g.nBlocks();
32
33 obj.order = order;
34 37
35 % Create the diffOps for each block 38 % Create the diffOps for each block
36 obj.diffOps = cell(1, nBlocks); 39 obj.diffOps = cell(1, nBlocks);
37 for i = 1:nBlocks 40 for i = 1:nBlocks
38 h = getHand(i); 41 h = getHand(i);
68 intf = g.connections{i,j}; 71 intf = g.connections{i,j};
69 if isempty(intf) 72 if isempty(intf)
70 continue 73 continue
71 end 74 end
72 75
73 76 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}, intfTypes{i,j});
74 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2});
75 D{i,i} = D{i,i} + ii; 77 D{i,i} = D{i,i} + ii;
76 D{i,j} = D{i,j} + ij; 78 D{i,j} = D{i,j} + ij;
77 79
78 [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); 80 [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}, intfTypes{i,j});
79 D{j,j} = D{j,j} + jj; 81 D{j,j} = D{j,j} + jj;
80 D{j,i} = D{j,i} + ji; 82 D{j,i} = D{j,i} + ji;
81 end 83 end
82 end 84 end
83 obj.D = blockmatrix.toMatrix(D); 85 obj.D = blockmatrix.toMatrix(D);
125 ops = sparse2cell(op, obj.NNN); 127 ops = sparse2cell(op, obj.NNN);
126 end 128 end
127 129
128 % Get a boundary operator specified by opName for the given boundary/BoundaryGroup 130 % Get a boundary operator specified by opName for the given boundary/BoundaryGroup
129 function op = getBoundaryOperator(obj, opName, boundary) 131 function op = getBoundaryOperator(obj, opName, boundary)
132
130 switch class(boundary) 133 switch class(boundary)
131 case 'cell' 134 case 'cell'
132 localOpName = [opName '_' boundary{2}];
133 blockId = boundary{1}; 135 blockId = boundary{1};
134 localOp = obj.diffOps{blockId}.(localOpName); 136 localOp = obj.diffOps{blockId}.getBoundaryOperator(opName, boundary{2});
135 137
136 div = {obj.blockmatrixDiv{1}, size(localOp,2)}; 138 div = {obj.blockmatrixDiv{1}, size(localOp,2)};
137 blockOp = blockmatrix.zero(div); 139 blockOp = blockmatrix.zero(div);
138 blockOp{blockId,1} = localOp; 140 blockOp{blockId,1} = localOp;
139 op = blockmatrix.toMatrix(blockOp); 141 op = blockmatrix.toMatrix(blockOp);
147 error('Unknown boundary indentifier') 149 error('Unknown boundary indentifier')
148 end 150 end
149 end 151 end
150 152
151 function op = getBoundaryQuadrature(obj, boundary) 153 function op = getBoundaryQuadrature(obj, boundary)
152 opName = 'H';
153 switch class(boundary) 154 switch class(boundary)
154 case 'cell' 155 case 'cell'
155 localOpName = [opName '_' boundary{2}];
156 blockId = boundary{1}; 156 blockId = boundary{1};
157 op = obj.diffOps{blockId}.(localOpName); 157 op = obj.diffOps{blockId}.getBoundaryQuadrature(boundary{2});
158
159 return 158 return
160 case 'multiblock.BoundaryGroup' 159 case 'multiblock.BoundaryGroup'
161 N = length(boundary); 160 N = length(boundary);
162 H_bm = cell(N,N); 161 H_bm = cell(N,N);
163 for i = 1:N 162 for i = 1:N
199 198
200 % Get the closure and penaly matrices 199 % Get the closure and penaly matrices
201 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); 200 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type);
202 201
203 % Expand to matrix for full domain. 202 % Expand to matrix for full domain.
204 div = obj.blockmatrixDiv; 203 closure = multiblock.local2globalClosure(blockClosure, obj.blockmatrixDiv, I);
205 if ~iscell(blockClosure) 204 penalty = multiblock.local2globalPenalty(blockPenalty, obj.blockmatrixDiv, I);
206 temp = blockmatrix.zero(div);
207 temp{I,I} = blockClosure;
208 closure = blockmatrix.toMatrix(temp);
209 else
210 for i = 1:length(blockClosure)
211 temp = blockmatrix.zero(div);
212 temp{I,I} = blockClosure{i};
213 closure{i} = blockmatrix.toMatrix(temp);
214 end
215 end
216
217 if ~iscell(blockPenalty)
218 div{2} = size(blockPenalty, 2); % Penalty is a column vector
219 p = blockmatrix.zero(div);
220 p{I} = blockPenalty;
221 penalty = blockmatrix.toMatrix(p);
222 else
223 % TODO: used by beam equation, should be eliminated. SHould only set one BC per call
224 for i = 1:length(blockPenalty)
225 div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector
226 p = blockmatrix.zero(div);
227 p{I} = blockPenalty{i};
228 penalty{i} = blockmatrix.toMatrix(p);
229 end
230 end
231 end 205 end
232 206
233 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 207 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
234 error('not implemented') 208 error('not implemented')
235 end 209 end