Mercurial > repos > public > sbplib

comparison +multiblock/DiffOp.m @ 781:69ab0e69f972 feature/interpolation

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/grids
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 24 Jul 2018 20:14:29 -0700
parents e05465aa2e25
children 5cf9fdf4c98f 517d6019d63c
comparison
equal deleted inserted replaced
751:005a8d071da3 781:69ab0e69f972
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, g, order, doParam)
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:
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, g, doParam);
30 30
31 nBlocks = grid.nBlocks(); 31 nBlocks = g.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);
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(g.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);
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 obj.blockmatrixDiv = {g.Ns, g.Ns};
57 D = blockmatrix.zero(obj.blockmatrixDiv); 57 D = blockmatrix.zero(obj.blockmatrixDiv);
58 for i = 1:nBlocks 58 for i = 1:nBlocks
59 D{i,i} = obj.diffOps{i}.D; 59 D{i,i} = obj.diffOps{i}.D;
60 end 60 end
61 61
62 for i = 1:nBlocks 62 for i = 1:nBlocks
63 for j = 1:nBlocks 63 for j = 1:nBlocks
64 intf = grid.connections{i,j}; 64 intf = g.connections{i,j};
65 if isempty(intf) 65 if isempty(intf)
66 continue 66 continue
67 end 67 end
68 68
69 69
75 D{j,j} = D{j,j} + jj; 75 D{j,j} = D{j,j} + jj;
76 D{j,i} = D{j,i} + ji; 76 D{j,i} = D{j,i} + ji;
77 end 77 end
78 end 78 end
79 obj.D = blockmatrix.toMatrix(D); 79 obj.D = blockmatrix.toMatrix(D);
80 80 obj.grid = g;
81 81
82 function [getHand, getParam] = parseInput(doHand, grid, doParam) 82
83 if ~isa(grid, 'multiblock.Grid') 83 function [getHand, getParam] = parseInput(doHand, g, doParam)
84 if ~isa(g, 'multiblock.Grid')
84 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); 85 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
85 end 86 end
86 87
87 if iscell(doHand) && length(doHand) == grid.nBlocks() 88 if iscell(doHand) && length(doHand) == g.nBlocks()
88 getHand = @(i)doHand{i}; 89 getHand = @(i)doHand{i};
89 elseif isa(doHand, 'function_handle') 90 elseif isa(doHand, 'function_handle')
90 getHand = @(i)doHand; 91 getHand = @(i)doHand;
91 else 92 else
92 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); 93 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks');
102 return 103 return
103 end 104 end
104 105
105 % doParam is a non-empty cell-array 106 % doParam is a non-empty cell-array
106 107
107 if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) 108 if length(doParam) == g.nBlocks() && all(cellfun(@iscell, doParam))
108 % doParam is a cell-array of cell-arrays 109 % doParam is a cell-array of cell-arrays
109 getParam = @(i)doParam{i}; 110 getParam = @(i)doParam{i};
110 return 111 return
111 end 112 end
112 113
114 end 115 end
115 end 116 end
116 117
117 function ops = splitOp(obj, op) 118 function ops = splitOp(obj, op)
118 % Splits a matrix operator into a cell-matrix of matrix operators for 119 % Splits a matrix operator into a cell-matrix of matrix operators for
119 % each grid. 120 % each g.
120 ops = sparse2cell(op, obj.NNN); 121 ops = sparse2cell(op, obj.NNN);
121 end 122 end
122 123
123 % Get a boundary operator specified by opName for the given boundary/BoundaryGroup 124 % Get a boundary operator specified by opName for the given boundary/BoundaryGroup
124 function op = getBoundaryOperator(obj, opName, boundary) 125 function op = getBoundaryOperator(obj, opName, boundary)
136 case 'multiblock.BoundaryGroup' 137 case 'multiblock.BoundaryGroup'
137 op = sparse(size(obj.D,1),0); 138 op = sparse(size(obj.D,1),0);
138 for i = 1:length(boundary) 139 for i = 1:length(boundary)
139 op = [op, obj.getBoundaryOperator(opName, boundary{i})]; 140 op = [op, obj.getBoundaryOperator(opName, boundary{i})];
140 end 141 end
142 otherwise
143 error('Unknown boundary indentifier')
144 end
145 end
146
147 function op = getBoundaryQuadrature(obj, boundary)
148 opName = 'H';
149 switch class(boundary)
150 case 'cell'
151 localOpName = [opName '_' boundary{2}];
152 blockId = boundary{1};
153 op = obj.diffOps{blockId}.(localOpName);
154
155 return
156 case 'multiblock.BoundaryGroup'
157 N = length(boundary);
158 H_bm = cell(N,N);
159 for i = 1:N
160 H_bm{i,i} = obj.getBoundaryQuadrature(boundary{i});
161 end
162 op = blockmatrix.toMatrix(H_bm);
141 otherwise 163 otherwise
142 error('Unknown boundary indentifier') 164 error('Unknown boundary indentifier')
143 end 165 end
144 end 166 end
145 167
192 div{2} = size(blockPenalty, 2); % Penalty is a column vector 214 div{2} = size(blockPenalty, 2); % Penalty is a column vector
193 p = blockmatrix.zero(div); 215 p = blockmatrix.zero(div);
194 p{I} = blockPenalty; 216 p{I} = blockPenalty;
195 penalty = blockmatrix.toMatrix(p); 217 penalty = blockmatrix.toMatrix(p);
196 else 218 else
219 % TODO: used by beam equation, should be eliminated. SHould only set one BC per call
197 for i = 1:length(blockPenalty) 220 for i = 1:length(blockPenalty)
198 div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector 221 div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector
199 p = blockmatrix.zero(div); 222 p = blockmatrix.zero(div);
200 p{I} = blockPenalty{i}; 223 p{I} = blockPenalty{i};
201 penalty{i} = blockmatrix.toMatrix(p); 224 penalty{i} = blockmatrix.toMatrix(p);