Mercurial > repos > public > sbplib

comparison +scheme/ViscoElastic2d.m @ 1312:6d64b57caf46 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add viscoelastic regular interfaces
author Martin Almquist <malmquist@stanford.edu>
date Tue, 21 Jul 2020 21:07:37 -0700
parents b2e84fe336d8
children a4c8bf2371f3
comparison
equal deleted inserted replaced
1311:b2e84fe336d8 1312:6d64b57caf46
335 defaultType.type = 'standard'; 335 defaultType.type = 'standard';
336 default_struct('type', defaultType); 336 default_struct('type', defaultType);
337 337
338 switch type.type 338 switch type.type
339 case 'standard' 339 case 'standard'
340 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type); 340 [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type);
341 case 'frictionalFault' 341 case 'frictionalFault'
342 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type); 342 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
343 end 343 end
344
345 end
346
347 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
348
349 % u denotes the solution in the own domain
350 % v denotes the solution in the neighbour domain
351 u = obj;
352 v = neighbour_scheme;
353
354 dim = obj.dim;
355
356 n = u.getNormal(boundary);
357 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
358 e = u.getBoundaryOperatorForScalarField('e', boundary);
359
360 ev = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
361
362 H = u.H;
363 RHO = u.RHO;
364 ETA = u.ETA;
365 C = u.C;
366 Eu = u.Eu;
367 eU = u.eU;
368 eGamma = u.eGamma;
369
370 CV = v.C;
371 Ev = v.Eu;
372 eV = v.eU;
373 eGammaV = v.eGamma;
374 nV = v.getNormal(neighbour_boundary);
375
376
377 % Get elastic closure and penalty
378 [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
379 closure = Eu*closure*Eu';
380 penalty = Eu*penalty*Ev';
381
382 % Add viscous part of traction coupling
383 for i = 1:dim
384 for j = 1:dim
385 for k = 1:dim
386 for l = 1:dim
387 closure = closure + 1/2*eU{j}*( (RHO*H)\(C{i,j,k,l}*e*H_gamma*n{i}*e'*eGamma{k,l}') );
388 penalty = penalty + 1/2*eU{j}*( (RHO*H)\(e*H_gamma*nV{i}*(ev'*CV{i,j,k,l}*ev)*ev'*eGammaV{k,l}') );
389 end
390 end
391 end
392 end
393
394 % Add penalty to strain rate eq
395 for i = 1:dim
396 for j = 1:dim
397 for k = 1:dim
398 for l = 1:dim
399 closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}*e'*eU{l}') );
400 penalty = penalty + 1/2*eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}*ev'*eV{l}') );
401 end
402 end
403 end
404 end
405
344 406
345 end 407 end
346 408
347 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type) 409 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
348 tuning = type.tuning; 410 tuning = type.tuning;