Mercurial > repos > public > sbplib

comparison +scheme/elasticDilationVariable.m @ 675:90bf651abc7c feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add Dirichlet BC to dilation. Stable with variable coeff and conv study yields p+1/2 which probably is ok.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 18 Jan 2018 13:36:56 -0800
parents dd84b8862aa8
children 9926efb39330
comparison
equal deleted inserted replaced
674:dd84b8862aa8 675:90bf651abc7c
170 obj.D = D; 170 obj.D = D;
171 %=========================================% 171 %=========================================%
172 172
173 % Divergence operator for BC 173 % Divergence operator for BC
174 Div = cell(dim,1); 174 Div = cell(dim,1);
175 % Loop over boundaries
175 for i = 1:dim 176 for i = 1:dim
176 Div{i} = sparse(m_tot,dim*m_tot); 177 Div{i} = sparse(m_tot,dim*m_tot);
178 % Loop over components
177 for j = 1:dim 179 for j = 1:dim
178 Div{i} = Div{i} + d(i,j)*(obj.e_l{i}*obj.d1_l{i}' + obj.e_r{i}*obj.d1_r{i}')*E{j}' ... 180 Div{i} = Div{i} + d(i,j)*(obj.e_l{i}*obj.d1_l{i}' + obj.e_r{i}*obj.d1_r{i}')*E{j}' ...
179 + db(i,j)*obj.D1{j}*E{j}'; 181 + db(i,j)*obj.D1{j}*E{j}';
180 end 182 end
181 end 183 end
232 h = obj.h; 234 h = obj.h;
233 235
234 switch type 236 switch type
235 % Dirichlet boundary condition 237 % Dirichlet boundary condition
236 case {'D','d','dirichlet'} 238 case {'D','d','dirichlet'}
237 error('Not implemented');
238 tuning = 1.2; 239 tuning = 1.2;
239 phi = obj.phi; 240 phi = obj.phi;
240 241
241 closures = cell(obj.dim,1); 242 % Initialize with zeros
242 penalties = cell(obj.dim,obj.dim); 243 m_tot = obj.grid.N();
243 % Loop over components 244 closure = sparse(m_tot*obj.dim, m_tot*obj.dim);
245 penalty = 0*E{j}*e{j};
246
247 % Loop over components to put penalties on
244 for i = 1:obj.dim 248 for i = 1:obj.dim
245 H_gamma_i = obj.H_boundary{i}; 249 sigma_i = tuning * obj.dim/(H11*h(j)) +...
246 sigma_ij = tuning*delta(i,j)*2/(gamma*h(j)) +... 250 tuning * 1/(H11*h(j)*phi);
247 tuning*delta_b(i,j)*(2/(H11*h(j)) + 2/(H11*h(j)*phi)); 251
248 252 closure = closure + E{i}*inv(RHO)*nj*Hi*...
249 ci = E{i}*inv(RHO)*nj*Hi*... 253 ( ...
250 ( (e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{i}' + ...
251 delta(i,j)*(e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{j}' ... 254 delta(i,j)*(e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{j}' ...
252 ) ... 255 + delta_b(i,j)*(e{j}*H_gamma*e{j}'*A*D1{i})'*E{j}' ...
253 - sigma_ij*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma*e{j}'*E{i}'; 256 ) ...
254 257 - delta(i,j)*sigma_i*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma*e{j}'*E{j}';
255 cj = E{j}*inv(RHO)*nj*Hi*... 258
256 ( delta_b(i,j)*(e{j}*H_gamma*e{j}'*A*D1{i})'*E{i}' ... 259 penalty = penalty - ...
257 ); 260 E{i}*inv(RHO)*nj*Hi*...
258 261 ( ...
259 if isempty(closures{i}) 262 delta(i,j)*(H_gamma*e{j}'*A*e{j}*d{j}')' ...
260 closures{i} = ci; 263 + delta_b(i,j)*(H_gamma*e{j}'*A*D1{i})' ...
261 else 264 ) ...
262 closures{i} = closures{i} + ci; 265 + delta(i,j)*sigma_i*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma;
263 end 266
264
265 if isempty(closures{j})
266 closures{j} = cj;
267 else
268 closures{j} = closures{j} + cj;
269 end
270
271 pii = - E{i}*inv(RHO)*nj*Hi*...
272 ( (H_gamma*e{j}'*A*e{j}*d{j}')' + ...
273 delta(i,j)*(H_gamma*e{j}'*A*e{j}*d{j}')' ...
274 ) ...
275 + sigma_ij*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma;
276
277 pji = - E{j}*inv(RHO)*nj*Hi*...
278 ( delta_b(i,j)*(H_gamma*e{j}'*A*D1{i})' );
279
280 % Dummies
281 pij = - 0*E{i}*e{j};
282 pjj = - 0*E{j}*e{j};
283
284 if isempty(penalties{i,i})
285 penalties{i,i} = pii;
286 else
287 penalties{i,i} = penalties{i,i} + pii;
288 end
289
290 if isempty(penalties{j,i})
291 penalties{j,i} = pji;
292 else
293 penalties{j,i} = penalties{j,i} + pji;
294 end
295
296 if isempty(penalties{i,j})
297 penalties{i,j} = pij;
298 else
299 penalties{i,j} = penalties{i,j} + pij;
300 end
301
302 if isempty(penalties{j,j})
303 penalties{j,j} = pjj;
304 else
305 penalties{j,j} = penalties{j,j} + pjj;
306 end
307 end 267 end
308 [rows, cols] = size(closures{1});
309 closure = sparse(rows, cols);
310 for i = 1:obj.dim
311 closure = closure + closures{i};
312 end
313 penalty = penalties;
314 268
315 % Free boundary condition 269 % Free boundary condition
316 case {'F','f','Free','free'} 270 case {'F','f','Free','free'}
317 closures = cell(obj.dim,1); 271 closures = cell(obj.dim,1);
318 penalties = cell(obj.dim,obj.dim); 272 penalties = cell(obj.dim,obj.dim);
319 273
320 % Divergence operator 274 % Divergence operator
321 Div = obj.Div{j}; 275 Div = obj.Div{j};
322 276
323 % Loop over components 277 closure = -nj*E{j}*inv(RHO)*Hi*e{j} ...
324 %for i = 1:obj.dim 278 *H_gamma*e{j}'*A*e{j}*e{j}'*Div;
325 closure = -nj*E{j}*inv(RHO)*Hi*e{j} ... 279 penalty = nj*E{j}*inv(RHO)*Hi*e{j} ...
326 *H_gamma*e{j}'*A*e{j}*e{j}'*Div; 280 *H_gamma*e{j}'*A*e{j};
327 penalty = nj*E{j}*inv(RHO)*Hi*e{j} ...
328 *H_gamma*e{j}'*A*e{j};
329 %end
330 % [rows, cols] = size(closures{1});
331 % closure = sparse(rows, cols);
332 % for i = 1:obj.dim
333 % closure = closure + closures{i};
334 % for j = 1:obj.dim
335 % if i~=j
336 % [rows cols] = size(penalties{j,j});
337 % penalties{i,j} = sparse(rows,cols);
338 % end
339 % end
340 % end
341 % penalty = penalties;
342 281
343 282
344 % Unknown boundary condition 283 % Unknown boundary condition
345 otherwise 284 otherwise
346 error('No such boundary condition: type = %s',type); 285 error('No such boundary condition: type = %s',type);