Mercurial > repos > public > sbplib

comparison +scheme/ViscoElastic2d.m @ 1310:eb015fe9605b feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
author Martin Almquist <malmquist@stanford.edu>
date Mon, 20 Jul 2020 22:06:33 -0700
parents f59b849df30f
children b2e84fe336d8
comparison
equal deleted inserted replaced
1309:f59b849df30f 1310:eb015fe9605b
290 case 2 290 case 2
291 dir = {0, 1}; 291 dir = {0, 1};
292 end 292 end
293 293
294 % Add viscous part of closure 294 % Add viscous part of closure
295 for j = 1:dim 295 for i = 1:dim
296 for i = 1:dim 296 for j = 1:dim
297 for k = 1:dim 297 for k = 1:dim
298 for l = 1:dim 298 for l = 1:dim
299 closure = closure + eU{j}*( (RHO*H)\(C{i,j,k,l}*e*dir{j}^2*H_gamma*n{i}*e'*eGamma{k,l}') ); 299 for m = 1:dim
300 closure = closure + eU{m}*( (RHO*H)\(C{i,j,k,l}*e*dir{m}*dir{j}*H_gamma*n{i}*e'*eGamma{k,l}') );
301 end
300 end 302 end
301 end 303 end
302 end 304 end
303 end 305 end
304 306
307 % Get elastic closure and penalty 309 % Get elastic closure and penalty
308 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning); 310 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning);
309 closure = Eu*closure*Eu'; 311 closure = Eu*closure*Eu';
310 penalty = Eu*penalty; 312 penalty = Eu*penalty;
311 313
314 switch component
315 case 't'
316 dir = obj.getTangent(boundary);
317 case 'n'
318 dir = obj.getNormal(boundary);
319 case 1
320 dir = {1, 0};
321 case 2
322 dir = {0, 1};
323 end
324
312 % Add penalty to strain rate eq 325 % Add penalty to strain rate eq
313 l = component;
314 for i = 1:dim 326 for i = 1:dim
315 for j = 1:dim 327 for j = 1:dim
316 for k = 1:dim 328 for k = 1:dim
317 closure = closure - eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}*e'*eU{l}') ); 329 for l = 1:dim
318 penalty = penalty + eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}) ); 330 for m = 1:dim
319 end 331 closure = closure - eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*dir{l}*dir{m}*n{k}*e'*eU{m}') );
320 end 332 end
321 end 333 penalty = penalty + eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*dir{l}*n{k}) );
322 334 end
323 end 335 end
324 336 end
325 end 337 end
326 338
327 function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning) 339 end
328 disp('WARNING: displacementBCNormalTangential is only guaranteed to work for displacement BC on one component and traction on the other'); 340
329 u = obj;
330
331 component = bc{1};
332 type = bc{2};
333
334 switch component
335 case 'n'
336 en = u.getBoundaryOperator('en', boundary);
337 tau_n = u.getBoundaryOperator('tau_n', boundary);
338 N = u.getNormal(boundary);
339 case 't'
340 en = u.getBoundaryOperator('et', boundary);
341 tau_n = u.getBoundaryOperator('tau_t', boundary);
342 N = u.getTangent(boundary);
343 end
344
345 % Operators
346 e = u.getBoundaryOperatorForScalarField('e', boundary);
347 h11 = u.getBorrowing(boundary);
348 n = u.getNormal(boundary);
349
350 C = u.C;
351 Ji = u.Ji;
352 s = spdiag(u.(['s_' boundary]));
353 m_tot = u.grid.N();
354
355 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
356 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
357
358 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
359 dim = u.dim;
360
361 % Preallocate
362 [~, m_int] = size(H_gamma);
363 closure = sparse(dim*m_tot, dim*m_tot);
364 penalty = sparse(dim*m_tot, m_int);
365
366 % Term 1: The symmetric term
367 Z = sparse(m_int, m_int);
368 for i = 1:dim
369 for j = 1:dim
370 for l = 1:dim
371 for k = 1:dim
372 Z = Z + n{i}*N{j}*e'*Ji*C{i,j,k,l}*e*n{k}*N{l};
373 end
374 end
375 end
376 end
377
378 Z = -tuning*dim*1/h11*s*Z;
379 closure = closure + en*H_gamma*Z*en';
380 penalty = penalty - en*H_gamma*Z;
381
382 % Term 2: The symmetrizing term
383 closure = closure + tau_n*H_gamma*en';
384 penalty = penalty - tau_n*H_gamma;
385
386 % Multiply all normal component terms by inverse of density x quadraure
387 closure = RHOi*Hi*closure;
388 penalty = RHOi*Hi*penalty;
389 end 341 end
390 342
391 % type Struct that specifies the interface coupling. 343 % type Struct that specifies the interface coupling.
392 % Fields: 344 % Fields:
393 % -- tuning: penalty strength, defaults to 1.0 345 % -- tuning: penalty strength, defaults to 1.0