Mercurial > repos > public > sbplib

comparison +multiblock/DiffOp.m @ 211:3c4ffbfbfb84 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged feature/grid into feature/beams
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 16 Jun 2016 10:56:47 +0200
parents 39b7dcb2c724
children 8b10476b9bb7
comparison
equal deleted inserted replaced
200:ef41fde95ac4 211:3c4ffbfbfb84
3 grid 3 grid
4 order 4 order
5 diffOps 5 diffOps
6 D 6 D
7 H 7 H
8
9 blockmatrixDiv
8 end 10 end
9 11
10 methods 12 methods
11 function obj = DiffOp(doHand, grid, order, doParam) 13 function obj = DiffOp(doHand, grid, order, doParam)
12 % 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
33 % Create the diffOps for each block 35 % Create the diffOps for each block
34 obj.diffOps = cell(1, nBlocks); 36 obj.diffOps = cell(1, nBlocks);
35 for i = 1:nBlocks 37 for i = 1:nBlocks
36 h = getHand(i); 38 h = getHand(i);
37 p = getParam(i); 39 p = getParam(i);
38 obj.diffOps{i} = h(grid.grid{i}, order, p{:}); 40 obj.diffOps{i} = h(grid.grids{i}, order, p{:});
39 end 41 end
40 42
41 43
42 % Build the norm matrix 44 % Build the norm matrix
43 H = cell(nBlocks, nBlocks); 45 H = cell(nBlocks, nBlocks);
44 for i = 1:nBlocks 46 for i = 1:nBlocks
45 H{i,i} = obj.diffOps{i}.H; 47 H{i,i} = obj.diffOps{i}.H;
46 end 48 end
47 obj.H = cell2sparse(H); 49 obj.H = blockmatrix.toMatrix(H);
48 50
49 51
50 % Build the differentiation matrix 52 % Build the differentiation matrix
51 D = cell(nBlocks, nBlocks); 53 D = cell(nBlocks, nBlocks);
52 for i = 1:nBlocks 54 for i = 1:nBlocks
58 intf = grid.connections{i,j}; 60 intf = grid.connections{i,j};
59 if isempty(intf) 61 if isempty(intf)
60 continue 62 continue
61 end 63 end
62 64
63 [ii, ij] = obj.diffOps{i}.inteface_coupling(intf{1}, obj.diffOps{j}, intf{2}); 65 [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2});
64 D{i,i} = D{i,i} + ii; 66 D{i,i} = D{i,i} + ii;
65 D{i,j} = D{i,j} + ij; 67 D{i,j} = D{i,j} + ij;
66 68
67 [jj, ji] = obj.diffOps{j}.inteface_coupling(intf{2}, obj.diffOps{i}, intf{1}); 69 [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1});
68 D{j,j} = D{j,j} + jj; 70 D{j,j} = D{j,j} + jj;
69 D{j,i} = D{j,i} + ji; 71 D{j,i} = D{j,i} + ji;
70 end 72 end
71 end 73 end
72 obj.D = cell2sparse(D); 74 obj.D = blockmatrix.toMatrix(D);
73 75
74 % end 76 obj.blockmatrixDiv = blockmatrix.getDivision(D);
75 77
76 function [getHand, getParam] = parseInput(doHand, grid, doParam) 78 function [getHand, getParam] = parseInput(doHand, grid, doParam)
77 if ~isa(grid, 'multiblock.Grid') 79 if ~isa(grid, 'multiblock.Grid')
78 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); 80 error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
79 end 81 end
100 % Splits a matrix operator into a cell-matrix of matrix operators for 102 % Splits a matrix operator into a cell-matrix of matrix operators for
101 % each grid. 103 % each grid.
102 ops = sparse2cell(op, obj.NNN); 104 ops = sparse2cell(op, obj.NNN);
103 end 105 end
104 106
105 function m = boundary_condition(obj,boundary,type,data) 107 % Creates the closere and penalty matrix for a given boundary condition,
108 % boundary -- the name of the boundary on the form [id][name] where
109 % id is the number of a block and name is the name of a
110 % boundary of that block example: 1s or 3w
111 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
112 I = boundary(1)
113 name = boundary(2:end);
114
115 % Get the closure and penaly matrices
116 [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type, data);
117
118 % Expand to matrix for full domain.
119 div = obj.blockmatrixDiv;
120 closure = blockmatrix.zero(div);
121 closure{I,I} = blockClosure;
122
123 div{2} = 1; % Penalty is a column vector
124 penalty = blockmatrix.zero(div);
125 penalty{I} = blockPenalty;
126
127 % Convert to matrix
128 closure = blockmatrix.toMatrix(closure);
129 penalty = blockmatrix.toMatrix(closure);
130 end
131
132 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
106 133
107 end 134 end
108 135
109 function m = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 136 % Size returns the number of degrees of freedom
110
111 end
112
113
114 function N = size(obj) 137 function N = size(obj)
115 N = 0; 138 N = 0;
116 for i = 1:length(diffOps) 139 for i = 1:length(obj.diffOps)
117 N = N + diffOps.size(); 140 N = N + obj.diffOps{i}.size();
118 end 141 end
119 end 142 end
120 end 143 end
121 end 144 end