Mercurial > repos > public > sbplib

comparison +scheme/Heat2dVariable.m @ 688:eb2f9233acc3 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add scheme Heat2dVariable
author Martin Almquist <malmquist@stanford.edu>
date Mon, 12 Feb 2018 13:39:39 -0800
parents
children 60eb7f46d8d9
comparison
equal deleted inserted replaced
687:e8fc3aa1faf6 688:eb2f9233acc3
1 classdef Heat2dVariable < scheme.Scheme
2
3 % Discretizes the Laplacian with variable coefficent,
4 % In the Heat equation way (i.e., the discretization matrix is not necessarily
5 % symmetric)
6 % u_t = div * (kappa * grad u )
7 % opSet should be cell array of opSets, one per dimension. This
8 % is useful if we have periodic BC in one direction.
9
10 properties
11 m % Number of points in each direction, possibly a vector
12 h % Grid spacing
13
14 grid
15 dim
16
17 order % Order of accuracy for the approximation
18
19 % Diagonal matrix for variable coefficients
20 KAPPA % Variable coefficient
21
22 D % Total operator
23 D1 % First derivatives
24
25 % Second derivatives
26 D2_kappa
27
28 H, Hi % Inner products
29 e_l, e_r
30 d1_l, d1_r % Normal derivatives at the boundary
31
32 H_boundary % Boundary inner products
33
34 end
35
36 methods
37
38 function obj = Heat2dVariable(g ,order, kappa_fun, opSet)
39 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
40 default_arg('kappa_fun', @(x,y) 0*x+1);
41 dim = 2;
42
43 assert(isa(g, 'grid.Cartesian'))
44
45 kappa = grid.evalOn(g, kappa_fun);
46 m = g.size();
47 m_tot = g.N();
48
49 h = g.scaling();
50 lim = g.lim;
51
52 % 1D operators
53 ops = cell(dim,1);
54 for i = 1:dim
55 ops{i} = opSet{i}(m(i), lim{i}, order);
56 end
57
58 I = cell(dim,1);
59 D1 = cell(dim,1);
60 D2 = cell(dim,1);
61 H = cell(dim,1);
62 Hi = cell(dim,1);
63 e_l = cell(dim,1);
64 e_r = cell(dim,1);
65 d1_l = cell(dim,1);
66 d1_r = cell(dim,1);
67
68 for i = 1:dim
69 I{i} = speye(m(i));
70 D1{i} = ops{i}.D1;
71 D2{i} = ops{i}.D2;
72 H{i} = ops{i}.H;
73 Hi{i} = ops{i}.HI;
74 e_l{i} = ops{i}.e_l;
75 e_r{i} = ops{i}.e_r;
76 d1_l{i} = ops{i}.d1_l;
77 d1_r{i} = ops{i}.d1_r;
78 end
79
80 %====== Assemble full operators ========
81 KAPPA = spdiag(kappa);
82 obj.KAPPA = KAPPA;
83
84 obj.D1 = cell(dim,1);
85 obj.D2_kappa = cell(dim,1);
86 obj.e_l = cell(dim,1);
87 obj.e_r = cell(dim,1);
88 obj.d1_l = cell(dim,1);
89 obj.d1_r = cell(dim,1);
90
91 % D1
92 obj.D1{1} = kron(D1{1},I{2});
93 obj.D1{2} = kron(I{1},D1{2});
94
95 % Boundary operators
96 obj.e_l{1} = kron(e_l{1},I{2});
97 obj.e_l{2} = kron(I{1},e_l{2});
98 obj.e_r{1} = kron(e_r{1},I{2});
99 obj.e_r{2} = kron(I{1},e_r{2});
100
101 obj.d1_l{1} = kron(d1_l{1},I{2});
102 obj.d1_l{2} = kron(I{1},d1_l{2});
103 obj.d1_r{1} = kron(d1_r{1},I{2});
104 obj.d1_r{2} = kron(I{1},d1_r{2});
105
106 % D2
107 for i = 1:dim
108 obj.D2_kappa{i} = sparse(m_tot);
109 end
110 ind = grid.funcToMatrix(g, 1:m_tot);
111
112 for i = 1:m(2)
113 D_kappa = D2{1}(kappa(ind(:,i)));
114 p = ind(:,i);
115 obj.D2_kappa{1}(p,p) = D_kappa;
116 end
117
118 for i = 1:m(1)
119 D_kappa = D2{2}(kappa(ind(i,:)));
120 p = ind(i,:);
121 obj.D2_kappa{2}(p,p) = D_kappa;
122 end
123
124 % Quadratures
125 obj.H = kron(H{1},H{2});
126 obj.Hi = inv(obj.H);
127 obj.H_boundary = cell(dim,1);
128 obj.H_boundary{1} = H{2};
129 obj.H_boundary{2} = H{1};
130
131 % Differentiation matrix D (without SAT)
132 D2_kappa = obj.D2_kappa;
133 D1 = obj.D1;
134 D = sparse(m_tot,m_tot);
135 for i = 1:dim
136 D = D + D2_kappa{i};
137 end
138 obj.D = D;
139 %=========================================%
140
141 % Misc.
142 obj.m = m;
143 obj.h = h;
144 obj.order = order;
145 obj.grid = g;
146 obj.dim = dim;
147
148 end
149
150
151 % Closure functions return the operators applied to the own domain to close the boundary
152 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
153 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
154 % type is a string specifying the type of boundary condition.
155 % data is a function returning the data that should be applied at the boundary.
156 % neighbour_scheme is an instance of Scheme that should be interfaced to.
157 % neighbour_boundary is a string specifying which boundary to interface to.
158 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
159 default_arg('type','Neumann');
160 default_arg('parameter', []);
161
162 % j is the coordinate direction of the boundary
163 % nj: outward unit normal component.
164 % nj = -1 for west, south, bottom boundaries
165 % nj = 1 for east, north, top boundaries
166 [j, nj] = obj.get_boundary_number(boundary);
167 switch nj
168 case 1
169 e = obj.e_r;
170 d = obj.d1_r;
171 case -1
172 e = obj.e_l;
173 d = obj.d1_l;
174 end
175
176 Hi = obj.Hi;
177 H_gamma = obj.H_boundary{j};
178 KAPPA = obj.KAPPA;
179 kappa_gamma = e{j}'*KAPPA*e{j};
180
181 switch type
182
183 % Dirichlet boundary condition
184 case {'D','d','dirichlet','Dirichlet'}
185 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' );
186 penalty = nj*Hi*d{j}*kappa_gamma*H_gamma;
187
188 % Free boundary condition
189 case {'N','n','neumann','Neumann'}
190 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' );
191 penalty = nj*Hi*e{j}*kappa_gamma*H_gamma;
192
193 % Unknown boundary condition
194 otherwise
195 error('No such boundary condition: type = %s',type);
196 end
197 end
198
199 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
200 % u denotes the solution in the own domain
201 % v denotes the solution in the neighbour domain
202 error('Interface not implemented');
203 end
204
205 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
206 function [j, nj] = get_boundary_number(obj, boundary)
207
208 switch boundary
209 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
210 j = 1;
211 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
212 j = 2;
213 otherwise
214 error('No such boundary: boundary = %s',boundary);
215 end
216
217 switch boundary
218 case {'w','W','west','West','s','S','south','South'}
219 nj = -1;
220 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
221 nj = 1;
222 end
223 end
224
225 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
226 function [return_op] = get_boundary_operator(obj, op, boundary)
227
228 switch boundary
229 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
230 j = 1;
231 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
232 j = 2;
233 otherwise
234 error('No such boundary: boundary = %s',boundary);
235 end
236
237 switch op
238 case 'e'
239 switch boundary
240 case {'w','W','west','West','s','S','south','South'}
241 return_op = obj.e_l{j};
242 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
243 return_op = obj.e_r{j};
244 end
245 case 'd'
246 switch boundary
247 case {'w','W','west','West','s','S','south','South'}
248 return_op = obj.d_l{j};
249 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
250 return_op = obj.d_r{j};
251 end
252 otherwise
253 error(['No such operator: operatr = ' op]);
254 end
255
256 end
257
258 function N = size(obj)
259 N = prod(obj.m);
260 end
261 end
262 end