Mercurial > repos > public > sbplib

comparison +multiblock/DiffOp.m @ 820:501750fbbfdb

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 Fri, 07 Sep 2018 14:40:58 +0200
parents 5cf9fdf4c98f
children 57760d7088ad 7d4f57725192
comparison
equal deleted inserted replaced
819:fdf0ef9150f4 820:501750fbbfdb
1 classdef DiffOp < scheme.Scheme
2 properties
3 grid
4 order
5 diffOps
6 D
7 H
8
9 blockmatrixDiv
10 end
11
12 methods
13 function obj = DiffOp(doHand, g, order, doParam)
14 % doHand -- may either be a function handle or a cell array of
15 % function handles for each grid. The function handle(s)
16 % should be on the form do = doHand(grid, order, ...)
17 % Additional parameters for each doHand may be provided in
18 % the doParam input.
19 % g -- a multiblock grid
20 % order -- integer specifying the order of accuracy
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
23 % to the number of blocks then each element is sent to the
24 % corresponding function handle as extra parameters:
25 % doHand(..., doParam{i}{:}) Otherwise doParam is sent as
26 % extra parameters to all doHand: doHand(..., doParam{:})
27 default_arg('doParam', [])
28
29 [getHand, getParam] = parseInput(doHand, g, doParam);
30
31 nBlocks = g.nBlocks();
32
33 obj.order = order;
34
35 % Create the diffOps for each block
36 obj.diffOps = cell(1, nBlocks);
37 for i = 1:nBlocks
38 h = getHand(i);
39 p = getParam(i);
40 if ~iscell(p)
41 p = {p};
42 end
43 obj.diffOps{i} = h(g.grids{i}, order, p{:});
44 end
45
46
47 % Build the norm matrix
48 H = cell(nBlocks, nBlocks);
49 for i = 1:nBlocks
50 H{i,i} = obj.diffOps{i}.H;
51 end
52 obj.H = blockmatrix.toMatrix(H);
53
54
55 % Build the differentiation matrix
56 Ns = zeros(nBlocks,1);
57 for i = 1:nBlocks
58 Ns(i) = length(obj.diffOps{i}.D);
59 end
60 obj.blockmatrixDiv = {Ns, Ns};
61 D = blockmatrix.zero(obj.blockmatrixDiv);
62 for i = 1:nBlocks
63 D{i,i} = obj.diffOps{i}.D;
64 end
65
66 for i = 1:nBlocks
67 for j = 1:nBlocks
68 intf = g.connections{i,j};
69 if isempty(intf)
70 continue
71 end
72
73
74 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2});
75 D{i,i} = D{i,i} + ii;
76 D{i,j} = D{i,j} + ij;
77
78 [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1});
79 D{j,j} = D{j,j} + jj;
80 D{j,i} = D{j,i} + ji;
81 end
82 end
83 obj.D = blockmatrix.toMatrix(D);
84 obj.grid = g;
85
86
87 function [getHand, getParam] = parseInput(doHand, g, doParam)
88 if ~isa(g, 'multiblock.Grid')
89 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
90 end
91
92 if iscell(doHand) && length(doHand) == g.nBlocks()
93 getHand = @(i)doHand{i};
94 elseif isa(doHand, 'function_handle')
95 getHand = @(i)doHand;
96 else
97 error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks');
98 end
99
100 if isempty(doParam)
101 getParam = @(i){};
102 return
103 end
104
105 if ~iscell(doParam)
106 getParam = @(i)doParam;
107 return
108 end
109
110 % doParam is a non-empty cell-array
111
112 if length(doParam) == g.nBlocks() && all(cellfun(@iscell, doParam))
113 % doParam is a cell-array of cell-arrays
114 getParam = @(i)doParam{i};
115 return
116 end
117
118 getParam = @(i)doParam;
119 end
120 end
121
122 function ops = splitOp(obj, op)
123 % Splits a matrix operator into a cell-matrix of matrix operators for
124 % each grid.
125 ops = sparse2cell(op, obj.NNN);
126 end
127
128 % Get a boundary operator specified by opName for the given boundary/BoundaryGroup
129 function op = getBoundaryOperator(obj, opName, boundary)
130 switch class(boundary)
131 case 'cell'
132 localOpName = [opName '_' boundary{2}];
133 blockId = boundary{1};
134 localOp = obj.diffOps{blockId}.(localOpName);
135
136 div = {obj.blockmatrixDiv{1}, size(localOp,2)};
137 blockOp = blockmatrix.zero(div);
138 blockOp{blockId,1} = localOp;
139 op = blockmatrix.toMatrix(blockOp);
140 return
141 case 'multiblock.BoundaryGroup'
142 op = sparse(size(obj.D,1),0);
143 for i = 1:length(boundary)
144 op = [op, obj.getBoundaryOperator(opName, boundary{i})];
145 end
146 otherwise
147 error('Unknown boundary indentifier')
148 end
149 end
150
151 function op = getBoundaryQuadrature(obj, boundary)
152 opName = 'H';
153 switch class(boundary)
154 case 'cell'
155 localOpName = [opName '_' boundary{2}];
156 blockId = boundary{1};
157 op = obj.diffOps{blockId}.(localOpName);
158
159 return
160 case 'multiblock.BoundaryGroup'
161 N = length(boundary);
162 H_bm = cell(N,N);
163 for i = 1:N
164 H_bm{i,i} = obj.getBoundaryQuadrature(boundary{i});
165 end
166 op = blockmatrix.toMatrix(H_bm);
167 otherwise
168 error('Unknown boundary indentifier')
169 end
170 end
171
172 % Creates the closure and penalty matrix for a given boundary condition,
173 % boundary -- the name of the boundary on the form {id,name} where
174 % id is the number of a block and name is the name of a
175 % boundary of that block example: {1,'s'} or {3,'w'}. It
176 % can also be a boundary group
177 function [closure, penalty] = boundary_condition(obj, boundary, type)
178 switch class(boundary)
179 case 'cell'
180 [closure, penalty] = obj.singleBoundaryCondition(boundary, type);
181 case 'multiblock.BoundaryGroup'
182 [n,m] = size(obj.D);
183 closure = sparse(n,m);
184 penalty = sparse(n,0);
185 for i = 1:length(boundary)
186 [closurePart, penaltyPart] = obj.boundary_condition(boundary{i}, type);
187 closure = closure + closurePart;
188 penalty = [penalty, penaltyPart];
189 end
190 otherwise
191 error('Unknown boundary indentifier')
192 end
193
194 end
195
196 function [closure, penalty] = singleBoundaryCondition(obj, boundary, type)
197 I = boundary{1};
198 name = boundary{2};
199
200 % Get the closure and penaly matrices
201 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type);
202
203 % Expand to matrix for full domain.
204 div = obj.blockmatrixDiv;
205 if ~iscell(blockClosure)
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
232
233 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
234 error('not implemented')
235 end
236
237 % Size returns the number of degrees of freedom
238 function N = size(obj)
239 N = 0;
240 for i = 1:length(obj.diffOps)
241 N = N + obj.diffOps{i}.size();
242 end
243 end
244 end
245 end