Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropic.m @ 1294:fe712a1fca3f feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add frictional fault interface in Elastic2dAnisotropicCurve
author Martin Almquist <malmquist@stanford.edu>
date Thu, 02 Jul 2020 15:00:21 -0700
parents a8e730db76e9
children cb053fabbedc
comparison
equal deleted inserted replaced
1293:3e2c1df740df 1294:fe712a1fca3f
33 33
34 % Traction operators 34 % Traction operators
35 tau_w, tau_e, tau_s, tau_n % Return vector field 35 tau_w, tau_e, tau_s, tau_n % Return vector field
36 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field 36 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
37 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field 37 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
38 tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field
38 39
39 % Inner products 40 % Inner products
40 H 41 H
41 42
42 % Boundary inner products (for scalar field) 43 % Boundary inner products (for scalar field)
317 obj.et_w = (obj.tangent_w{1}*obj.e1_w' + obj.tangent_w{2}*obj.e2_w')'; 318 obj.et_w = (obj.tangent_w{1}*obj.e1_w' + obj.tangent_w{2}*obj.e2_w')';
318 obj.et_e = (obj.tangent_e{1}*obj.e1_e' + obj.tangent_e{2}*obj.e2_e')'; 319 obj.et_e = (obj.tangent_e{1}*obj.e1_e' + obj.tangent_e{2}*obj.e2_e')';
319 obj.et_s = (obj.tangent_s{1}*obj.e1_s' + obj.tangent_s{2}*obj.e2_s')'; 320 obj.et_s = (obj.tangent_s{1}*obj.e1_s' + obj.tangent_s{2}*obj.e2_s')';
320 obj.et_n = (obj.tangent_n{1}*obj.e1_n' + obj.tangent_n{2}*obj.e2_n')'; 321 obj.et_n = (obj.tangent_n{1}*obj.e1_n' + obj.tangent_n{2}*obj.e2_n')';
321 322
323 % obj.tau_n_w = (obj.en_w'*obj.e_w*obj.tau_w')';
324 % obj.tau_n_e = (obj.en_e'*obj.e_e*obj.tau_e')';
325 % obj.tau_n_s = (obj.en_s'*obj.e_s*obj.tau_s')';
326 % obj.tau_n_n = (obj.en_n'*obj.e_n*obj.tau_n')';
327
328 obj.tau_n_w = (obj.n_w{1}*obj.tau1_w' + obj.n_w{2}*obj.tau2_w')';
329 obj.tau_n_e = (obj.n_e{1}*obj.tau1_e' + obj.n_e{2}*obj.tau2_e')';
330 obj.tau_n_s = (obj.n_s{1}*obj.tau1_s' + obj.n_s{2}*obj.tau2_s')';
331 obj.tau_n_n = (obj.n_n{1}*obj.tau1_n' + obj.n_n{2}*obj.tau2_n')';
332
322 % Stress operators 333 % Stress operators
323 sigma = cell(dim, dim); 334 sigma = cell(dim, dim);
324 D1 = {obj.Dx, obj.Dy}; 335 D1 = {obj.Dx, obj.Dy};
325 E = obj.E; 336 E = obj.E;
326 N = length(obj.RHO); 337 N = length(obj.RHO);
410 penalty = penalty*s; 421 penalty = penalty*s;
411 422
412 end 423 end
413 end 424 end
414 425
426 function [closure, penalty] = displacementBCNormalTangential(boundary, bc, tuning)
427 error('not implemented');
428 u = obj;
429
430 component = bc{1};
431 type = bc{2};
432
433 switch component
434 case 'n'
435 en_u = u.getBoundaryOperator('en', boundary);
436 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
437 case 't'
438 en_u = u.getBoundaryOperator('et', boundary);
439 tau_n_u = u.getBoundaryOperator('tau_t', boundary);
440 end
441
442 % Operators
443 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
444 h11_u = u.getBorrowing(boundary);
445 n_u = u.getNormal(boundary);
446
447 C_u = u.C;
448 Ji_u = u.Ji;
449 s_u = spdiag(u.(['s_' boundary]));
450 m_tot_u = u.grid.N();
451
452 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
453 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
454
455 % Shared operators
456 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
457 dim = u.dim;
458
459 % Preallocate
460 [~, m_int] = size(H_gamma);
461 closure = sparse(dim*m_tot_u, dim*m_tot_u);
462 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
463
464 % Continuity of normal displacement, term 1: The symmetric term
465 Z_u = sparse(m_int, m_int);
466 Z_v = sparse(m_int, m_int);
467 for i = 1:dim
468 for j = 1:dim
469 for l = 1:dim
470 for k = 1:dim
471 Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{l,i,k,j}*e_u*n_u{k}*n_u{l};
472 Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{l,i,k,j}*e_v*n_v{k}*n_v{l};
473 end
474 end
475 end
476 end
477
478 Z = -tuning*dim*( 1/(4*h11_u)*s_u*Z_u + 1/(4*h11_v)*s_v*Z_v );
479 closure = closure + en_u*H_gamma*Z*en_u';
480 penalty = penalty + en_u*H_gamma*Z*en_v';
481
482 % Continuity of normal displacement, term 2: The symmetrizing term
483 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
484 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
485
486 % Continuity of normal traction
487 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
488 penalty = penalty - 1/2*en_u*H_gamma*tau_n_v';
489
490 % Multiply all normal component terms by inverse of density x quadraure
491 closure = RHOi*Hi*closure;
492 penalty = RHOi*Hi*penalty;
493 end
494
415 % type Struct that specifies the interface coupling. 495 % type Struct that specifies the interface coupling.
416 % Fields: 496 % Fields:
417 % -- tuning: penalty strength, defaults to 1.0 497 % -- tuning: penalty strength, defaults to 1.0
418 % -- interpolation: type of interpolation, default 'none' 498 % -- interpolation: type of interpolation, default 'none'
419 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) 499 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
420 500
421 defaultType.tuning = 1.0; 501 defaultType.tuning = 1.0;
422 defaultType.interpolation = 'none'; 502 defaultType.interpolation = 'none';
503 defaultType.type = 'standard';
423 default_struct('type', defaultType); 504 default_struct('type', defaultType);
424 505
425 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type); 506 switch type.type
426 end 507 case 'standard'
508 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
509 case 'frictionalFault'
510 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
511 end
512
513 end
514
515 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
516 tuning = type.tuning;
517
518 % u denotes the solution in the own domain
519 % v denotes the solution in the neighbour domain
520
521 u = obj;
522 v = neighbour_scheme;
523
524 % Operators, u side
525 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
526 en_u = u.getBoundaryOperator('en', boundary);
527 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
528 h11_u = u.getBorrowing(boundary);
529 n_u = u.getNormal(boundary);
530
531 C_u = u.C;
532 Ji_u = u.Ji;
533 s_u = spdiag(u.(['s_' boundary]));
534 m_tot_u = u.grid.N();
535
536 % Operators, v side
537 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
538 en_v = v.getBoundaryOperator('en', neighbour_boundary);
539 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
540 h11_v = v.getBorrowing(neighbour_boundary);
541 n_v = v.getNormal(neighbour_boundary);
542
543 C_v = v.C;
544 Ji_v = v.Ji;
545 s_v = spdiag(v.(['s_' neighbour_boundary]));
546 m_tot_v = v.grid.N();
547
548 % Operators that are only required for own domain
549 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
550 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
551
552 % Shared operators
553 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
554 dim = u.dim;
555
556 % Preallocate
557 [~, m_int] = size(H_gamma);
558 closure = sparse(dim*m_tot_u, dim*m_tot_u);
559 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
560
561 % Continuity of normal displacement, term 1: The symmetric term
562 Z_u = sparse(m_int, m_int);
563 Z_v = sparse(m_int, m_int);
564 for i = 1:dim
565 for j = 1:dim
566 for l = 1:dim
567 for k = 1:dim
568 Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{l,i,k,j}*e_u*n_u{k}*n_u{l};
569 Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{l,i,k,j}*e_v*n_v{k}*n_v{l};
570 end
571 end
572 end
573 end
574
575 Z = -tuning*dim*( 1/(4*h11_u)*s_u*Z_u + 1/(4*h11_v)*s_v*Z_v );
576 closure = closure + en_u*H_gamma*Z*en_u';
577 penalty = penalty + en_u*H_gamma*Z*en_v';
578
579 % Continuity of normal displacement, term 2: The symmetrizing term
580 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
581 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
582
583 % Continuity of normal traction
584 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
585 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
586
587 % Multiply all normal component terms by inverse of density x quadraure
588 closure = RHOi*Hi*closure;
589 penalty = RHOi*Hi*penalty;
590
591 % ---- Tangential tractions are imposed just like traction BC ------
592 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
593
594 end
595
427 596
428 % Returns h11 for the boundary specified by the string boundary. 597 % Returns h11 for the boundary specified by the string boundary.
429 % op -- string 598 % op -- string
430 function h11 = getBorrowing(obj, boundary) 599 function h11 = getBorrowing(obj, boundary)
431 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 600 assertIsMember(boundary, {'w', 'e', 's', 'n'})
448 617
449 % Returns the boundary operator op for the boundary specified by the string boundary. 618 % Returns the boundary operator op for the boundary specified by the string boundary.
450 % op -- string 619 % op -- string
451 function o = getBoundaryOperator(obj, op, boundary) 620 function o = getBoundaryOperator(obj, op, boundary)
452 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 621 assertIsMember(boundary, {'w', 'e', 's', 'n'})
453 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et'}) 622 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n'})
454 623
455 o = obj.([op, '_', boundary]); 624 o = obj.([op, '_', boundary]);
456 625
457 end 626 end
458 627