Mercurial > repos > public > sbplib

comparison +scheme/bcSetup.m @ 890:c70131daaa6e feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/poroelastic.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 21 Nov 2018 18:29:29 -0800
parents 675e571e6f4a 76efb6a7b466
children 386ef449df51
comparison
equal deleted inserted replaced
885:18e10217dca9 890:c70131daaa6e
9 % 9 %
10 % The boundary data function can either be a function of time or a function of time and space coordinates. 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. 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. 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). 13 % For example in the 2D case: f(t,x,y).
14 function [closure, S] = bcSetup(diffOp, bc, S_sign) 14 function [closure, S] = bcSetup(diffOp, bcs, S_sign)
15 default_arg('S_sign', 1); 15 default_arg('S_sign', 1);
16 assertType(bc, 'cell'); 16 assertType(bcs, 'cell');
17 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');
18 18
19 verifyBcFormat(bcs, diffOp);
19 20
20 % Setup storage arrays 21 % Setup storage arrays
21 closure = spzeros(size(diffOp)); 22 closure = spzeros(size(diffOp));
22 gridDataPenalties = {}; 23 gridData = {};
23 gridDataFunctions = {}; 24 symbolicData = {};
24 symbolicDataPenalties = {};
25 symbolicDataFunctions = {};
26 symbolicDataCoords = {};
27 25
28 % Collect closures, penalties and data 26 % Collect closures, penalties and data
29 for i = 1:length(bc) 27 for i = 1:length(bcs)
30 assertType(bc{i}, 'struct'); 28 [localClosure, penalty] = diffOp.boundary_condition(bcs{i}.boundary, bcs{i}.type);
31 [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type);
32 closure = closure + localClosure; 29 closure = closure + localClosure;
33 30
34 if ~isfield(bc{i},'data') || isempty(bc{i}.data) 31 [ok, isSymbolic, data] = parseData(bcs{i}, penalty, diffOp.grid);
35 % Skip to next loop if there is no data 32
33 if ~ok
34 % There was no data
36 continue 35 continue
37 end 36 end
38 assertType(bc{i}.data, 'function_handle');
39 37
40 % Find dimension 38 if isSymbolic
41 dim = size(diffOp.grid.getBoundary(bc{i}.boundary), 2); 39 symbolicData{end+1} = data;
42
43 if nargin(bc{i}.data) == 1
44 % Grid data
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 40 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); 41 gridData{end+1} = data;
57 end 42 end
58 end 43 end
59 44
60 % Setup penalty function 45 % Setup penalty function
61 O = spzeros(size(diffOp),1); 46 O = spzeros(size(diffOp),1);
62 function v = S_fun(t) 47 function v = S_fun(t)
63 v = O; 48 v = O;
64 for i = 1:length(gridDataFunctions) 49 for i = 1:length(gridData)
65 v = v + gridDataPenalties{i}*gridDataFunctions{i}(t); 50 v = v + gridData{i}.penalty*gridData{i}.func(t);
66 end 51 end
67 52
68 for i = 1:length(symbolicDataFunctions) 53 for i = 1:length(symbolicData)
69 v = v + symbolicDataPenalties{i}*symbolicDataFunctions{i}(t, symbolicDataCoords{i}{:}); 54 v = v + symbolicData{i}.penalty*symbolicData{i}.func(t, symbolicData{i}.coords{:});
70 end 55 end
71 56
72 v = S_sign * v; 57 v = S_sign * v;
73 end 58 end
74 S = @S_fun; 59 S = @S_fun;
75 end 60 end
61
62 function verifyBcFormat(bcs, diffOp)
63 for i = 1:length(bcs)
64 assertType(bcs{i}, 'struct');
65 assertStructFields(bcs{i}, {'type', 'boundary'});
66
67 if ~isfield(bcs{i}, 'data') || isempty(bcs{i}.data)
68 continue
69 end
70
71 if ~isa(bcs{i}.data, 'function_handle')
72 error('bcs{%d}.data should be a function of time or a function of time and space',i);
73 end
74
75 b = diffOp.grid.getBoundary(bcs{i}.boundary);
76
77 dim = size(b,2);
78
79 if nargin(bcs{i}.data) == 1
80 % Grid data (only function of time)
81 assertSize(bcs{i}.data(0), 1, size(b));
82 elseif nargin(bcs{i}.data) ~= 1+dim
83 error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bcs{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i);
84 end
85 end
86 end
87
88 function [ok, isSymbolic, dataStruct] = parseData(bc, penalty, grid)
89 if ~isfield(bc,'data') || isempty(bc.data)
90 isSymbolic = [];
91 dataStruct = struct();
92 ok = false;
93 return
94 end
95 ok = true;
96
97 nArg = nargin(bc.data);
98
99 if nArg > 1
100 % Symbolic data
101 isSymbolic = true;
102 coord = grid.getBoundary(bc.boundary);
103 dataStruct.penalty = penalty;
104 dataStruct.func = bc.data;
105 dataStruct.coords = num2cell(coord, 1);
106 else
107 % Grid data
108 isSymbolic = false;
109 dataStruct.penalty = penalty;
110 dataStruct.func = bcs{i}.data;
111 end
112 end