Mercurial > repos > public > sbplib
comparison +multiblock/stitchSchemes.m @ 0:48b6fb693025
Initial commit.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Thu, 17 Sep 2015 10:12:50 +0200 |
| parents | |
| children | 5c569cbef49e |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:48b6fb693025 |
|---|---|
| 1 % Stitch schemes together given connection matrix and BC vector. | |
| 2 % schmHand - function_handle to a Scheme constructor | |
| 3 % order - order of accuracy | |
| 4 % schmParam - cell array of extra parameters sent to each Scheme stored as cell arrays | |
| 5 % blocks - block definitions, On whatever form the scheme expects. | |
| 6 % ms - grid points in each direction for each block. Ex {[10,10], [10, 20]} | |
| 7 % conn - connection matrix | |
| 8 % bound - boundary condition vector, array of structs with fields w,e,s,n | |
| 9 % each field with a parameter array that is sent to schm.boundary_condition | |
| 10 % | |
| 11 % Output parameters are cell arrays and cell matrices. | |
| 12 function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound) | |
| 13 | |
| 14 n_blocks = numel(blocks); | |
| 15 | |
| 16 % Creating Schemes | |
| 17 for i = 1:n_blocks | |
| 18 if ~iscell(schmParam{i}) | |
| 19 param = schmParam(i); | |
| 20 else | |
| 21 param = schmParam{i}; | |
| 22 end | |
| 23 | |
| 24 % class(schmParam) | |
| 25 % class(ms) | |
| 26 % class(blocks) | |
| 27 % class(schmParam{i}) | |
| 28 % class(ms) | |
| 29 | |
| 30 schms{i} = schmHand(ms{i},blocks{i},order,param{:}); | |
| 31 end | |
| 32 | |
| 33 | |
| 34 % Total norm | |
| 35 H = cell(n_blocks,n_blocks); | |
| 36 for i = 1:n_blocks | |
| 37 H{i,i} = schms{i}.H; | |
| 38 end | |
| 39 | |
| 40 | |
| 41 | |
| 42 | |
| 43 | |
| 44 %% Total system matrix | |
| 45 | |
| 46 % Differentiation terms | |
| 47 D = cell(n_blocks,n_blocks); | |
| 48 for i = 1:n_blocks | |
| 49 D{i,i} = schms{i}.D; | |
| 50 end | |
| 51 | |
| 52 % Boundary penalty terms | |
| 53 for i = 1:n_blocks | |
| 54 if ~isstruct(bound{i}) | |
| 55 continue | |
| 56 end | |
| 57 | |
| 58 fn = fieldnames(bound{i}); | |
| 59 for j = 1:length(fn); | |
| 60 bc = bound{i}.(fn{j}); | |
| 61 if isempty(bc) | |
| 62 continue | |
| 63 end | |
| 64 | |
| 65 t = schms{i}.boundary_condition(fn{j},bc{:}); | |
| 66 D{i,i} = D{i,i}+t; | |
| 67 end | |
| 68 end | |
| 69 | |
| 70 % Interface penalty terms | |
| 71 for i = 1:n_blocks | |
| 72 for j = 1:n_blocks | |
| 73 intf = conn{i,j}; | |
| 74 if isempty(intf) | |
| 75 continue | |
| 76 end | |
| 77 | |
| 78 [uu,uv,vv,vu] = noname.Scheme.interface_coupling(schms{i},intf{1},schms{j},intf{2}); | |
| 79 D{i,i} = D{i,i} + uu; | |
| 80 D{i,j} = uv; | |
| 81 D{j,j} = D{j,j} + vv; | |
| 82 D{j,i} = vu; | |
| 83 end | |
| 84 end | |
| 85 end |