Mercurial > repos > public > sbplib

comparison +scheme/bcSetup.m @ 749:1de60c4d462d feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Allow gird functions as data in bcSetup
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 04 Jun 2018 17:47:32 +0200
parents 87436a107d8a
children 2645188489f6 0776fa4754ff
comparison
equal deleted inserted replaced
748:ec0a594bea9f 749:1de60c4d462d
1 % function [closure, S] = bcSetup(diffOp, bc) 1 % function [closure, S] = bcSetup(diffOp, bc)
2 % Takes a diffOp and a cell array of boundary condition definitions. 2 % Takes a diffOp and a cell array of boundary condition definitions.
3 % Each bc is a struct with the fields 3 % Each bc is a struct with the fields
4 % * type -- Type of boundary condition 4 % * type -- Type of boundary condition
5 % * boundary -- Boundary identifier 5 % * boundary -- Boundary identifier
6 % * data -- A function_handle with time and space coordinates as a parameters, for example f(t,x,y) for a 2D problem 6 % * data -- A function_handle for a function which provides boundary data.(see below)
7 % Also takes S_sign which modifies the sign of S, [-1,1] 7 % Also takes S_sign which modifies the sign of S, [-1,1]
8 % Returns a closure matrix and a forcing function S 8 % Returns a closure matrix and a forcing function S.
9 %
10 % The boundary data function can either be a function of time or a function of time and space coordinates.
11 % In the case where it only depends on time it should return the data as grid function for the boundary.
12 % In the case where it also takes space coordinates the number of space coordinates should match the number of dimensions of the problem domain.
13 % For example in the 2D case: f(t,x,y).
9 function [closure, S] = bcSetup(diffOp, bc, S_sign) 14 function [closure, S] = bcSetup(diffOp, bc, S_sign)
10 default_arg('S_sign', 1); 15 default_arg('S_sign', 1);
11 assertType(bc, 'cell'); 16 assertType(bc, 'cell');
12 assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1'); 17 assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1');
13 18
14 19
20 % Setup storage arrays
15 closure = spzeros(size(diffOp)); 21 closure = spzeros(size(diffOp));
16 penalties = {}; 22 gridDataPenalties = {};
17 dataFunctions = {}; 23 gridDataFunctions = {};
18 dataParams = {}; 24 symbolicDataPenalties = {};
25 symbolicDataFunctions = {};
26 symbolicDataCoords = {};
19 27
28 % Collect closures, penalties and data
20 for i = 1:length(bc) 29 for i = 1:length(bc)
21 assertType(bc{i}, 'struct'); 30 assertType(bc{i}, 'struct');
22 [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type); 31 [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type);
23 closure = closure + localClosure; 32 closure = closure + localClosure;
24 33
25 if ~isfield(bc{i},'data') || isempty(bc{i}.data) 34 if ~isfield(bc{i},'data') || isempty(bc{i}.data)
35 % Skip to next loop if there is no data
26 continue 36 continue
27 end 37 end
28 assertType(bc{i}.data, 'function_handle'); 38 assertType(bc{i}.data, 'function_handle');
29 39
30 coord = diffOp.grid.getBoundary(bc{i}.boundary); 40 % Find dimension
31 assertNumberOfArguments(bc{i}.data, 1+size(coord,2)); 41 dim = size(diffOp.grid.getBoundary(bc{i}.boundary), 2);
32 42
33 penalties{end+1} = penalty; 43 if nargin(bc{i}.data) == 1
34 dataFunctions{end+1} = bc{i}.data; 44 % Grid data
35 dataParams{end+1} = num2cell(coord ,1); 45 boundarySize = [size(diffOp.grid.getBoundary(bc{i}.boundary),1),1];
46 assert_size(bc{i}.data(0), boundarySize); % Eval for t = 0 and make sure the function returns a grid vector of the correct size.
47 gridDataPenalties{end+1} = penalty;
48 gridDataFunctions{end+1} = bc{i}.data;
49 elseif nargin(bc{i}.data) == 1+dim
50 % Symbolic data
51 coord = diffOp.grid.getBoundary(bc{i}.boundary);
52 symbolicDataPenalties{end+1} = penalty;
53 symbolicDataFunctions{end+1} = bc{i}.data;
54 symbolicDataCoords{end+1} = num2cell(coord ,1);
55 else
56 error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bc{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i);
57 end
36 end 58 end
37 59
60 % Setup penalty function
38 O = spzeros(size(diffOp),1); 61 O = spzeros(size(diffOp),1);
39 function v = S_fun(t) 62 function v = S_fun(t)
40 v = O; 63 v = O;
41 for i = 1:length(dataFunctions) 64 for i = 1:length(gridDataFunctions)
42 v = v + penalties{i}*dataFunctions{i}(t, dataParams{i}{:}); 65 v = v + gridDataPenalties{i}*gridDataFunctions{i}(t);
66 end
67
68 for i = 1:length(symbolicDataFunctions)
69 v = v + symbolicDataPenalties{i}*symbolicDataFunctions{i}(t, symbolicDataCoords{i}{:});
43 end 70 end
44 71
45 v = S_sign * v; 72 v = S_sign * v;
46 end 73 end
47 S = @S_fun; 74 S = @S_fun;