Mercurial > repos > public > sbplib

comparison +multiblock/DiffOp.m @ 1331:60c875c18de3 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/poroelastic for Elastic schemes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 10 Mar 2022 16:54:26 +0100
parents 9c8ed00732fd bd8e607768ce
children
comparison
equal deleted inserted replaced
1330:855871e0b852 1331:60c875c18de3
127 ops = sparse2cell(op, obj.NNN); 127 ops = sparse2cell(op, obj.NNN);
128 end 128 end
129 129
130 % 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
131 function op = getBoundaryOperator(obj, opName, boundary) 131 function op = getBoundaryOperator(obj, opName, boundary)
132 blockmatrixDiv = obj.blockmatrixDiv{1};
132 133
133 switch class(boundary) 134 switch class(boundary)
134 case 'cell' 135 case 'cell'
135 blockId = boundary{1}; 136 blockId = boundary{1};
136 localOp = obj.diffOps{blockId}.getBoundaryOperator(opName, boundary{2}); 137 localOp = obj.diffOps{blockId}.getBoundaryOperator(opName, boundary{2});
137 138
138 div = {obj.blockmatrixDiv{1}, size(localOp,2)}; 139 div = {blockmatrixDiv, size(localOp,2)};
139 blockOp = blockmatrix.zero(div); 140 blockOp = blockmatrix.zero(div);
140 blockOp{blockId,1} = localOp; 141 blockOp{blockId,1} = localOp;
141 op = blockmatrix.toMatrix(blockOp); 142 op = blockmatrix.toMatrix(blockOp);
142 return 143 return
143 case 'multiblock.BoundaryGroup' 144 case 'multiblock.BoundaryGroup'
144 op = sparse(size(obj.D,1),0); 145 op = sparse(sum(blockmatrixDiv),0);
145 for i = 1:length(boundary) 146 for i = 1:length(boundary)
146 op = [op, obj.getBoundaryOperator(opName, boundary{i})]; 147 op = [op, obj.getBoundaryOperator(opName, boundary{i})];
147 end 148 end
148 otherwise 149 otherwise
149 error('Unknown boundary indentifier') 150 error('Unknown boundary indentifier')
206 207
207 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 208 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
208 error('not implemented') 209 error('not implemented')
209 end 210 end
210 211
212 function penalties = interfaceForcing(obj, boundary, neighbour_boundary, type)
213 default_arg('type', []);
214 switch class(boundary)
215 case 'cell'
216 penalties = obj.singleInterfaceForcing(boundary, neighbour_boundary, type);
217 case 'multiblock.BoundaryGroup'
218 [n,m] = size(obj.D);
219 penaltyZero = sparse(n,0);
220
221 for i = 1:length(boundary)
222 penaltyParts = obj.interfaceForcing(boundary{i}, neighbour_boundary{i}, type);
223
224 % Initalize if this is the first boundary
225 if i==1
226 penalties = cell(numel(penaltyParts), 1);
227 for j = 1:numel(penaltyParts)
228 penalties{j} = penaltyZero;
229 end
230 end
231
232 for j = 1:numel(penaltyParts)
233 penalties{j} = [penalties{j}, penaltyParts{j}];
234 end
235
236 end
237 otherwise
238 error('Unknown boundary indentifier')
239 end
240
241 end
242
243 function penalties = singleInterfaceForcing(obj, boundary, neighbour_boundary, type)
244 I = boundary{1};
245 b1 = boundary{2};
246
247 J = neighbour_boundary{1};
248 b2 = neighbour_boundary{2};
249
250 % Get local penalties
251 [~, ~, blockPenalties] = obj.diffOps{I}.interface(b1, obj.diffOps{J}, b2, type);
252 % [~, ~, blockPenalties2] = obj.diffOps{J}.interface(b2, b1, obj.diffOps{I});
253
254 % Expand to matrices for full domain.
255 n = numel(blockPenalties);
256 penalties = cell(n, 1);
257 for i = 1:n
258 penalties{i} = multiblock.local2globalPenalty(blockPenalties{i}, obj.blockmatrixDiv, I);
259 % + multiblock.local2globalPenalty(blockPenalties2{i}, obj.blockmatrixDiv, J);
260 end
261 end
262
211 % Size returns the number of degrees of freedom 263 % Size returns the number of degrees of freedom
212 function N = size(obj) 264 function N = size(obj)
213 N = 0; 265 N = 0;
214 for i = 1:length(obj.diffOps) 266 for i = 1:length(obj.diffOps)
215 N = N + obj.diffOps{i}.size(); 267 N = N + obj.diffOps{i}.size();