Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropic.m @ 1295:cb053fabbedc feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 02 Jul 2020 20:05:50 -0700
parents fe712a1fca3f
children a0d615bde7f8
comparison
equal deleted inserted replaced
1294:fe712a1fca3f 1295:cb053fabbedc
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 tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field
39 tau_t_w, tau_t_e, tau_t_s, tau_t_n % Return scalar field
39 40
40 % Inner products 41 % Inner products
41 H 42 H
42 43
43 % Boundary inner products (for scalar field) 44 % Boundary inner products (for scalar field)
318 obj.et_w = (obj.tangent_w{1}*obj.e1_w' + obj.tangent_w{2}*obj.e2_w')'; 319 obj.et_w = (obj.tangent_w{1}*obj.e1_w' + obj.tangent_w{2}*obj.e2_w')';
319 obj.et_e = (obj.tangent_e{1}*obj.e1_e' + obj.tangent_e{2}*obj.e2_e')'; 320 obj.et_e = (obj.tangent_e{1}*obj.e1_e' + obj.tangent_e{2}*obj.e2_e')';
320 obj.et_s = (obj.tangent_s{1}*obj.e1_s' + obj.tangent_s{2}*obj.e2_s')'; 321 obj.et_s = (obj.tangent_s{1}*obj.e1_s' + obj.tangent_s{2}*obj.e2_s')';
321 obj.et_n = (obj.tangent_n{1}*obj.e1_n' + obj.tangent_n{2}*obj.e2_n')'; 322 obj.et_n = (obj.tangent_n{1}*obj.e1_n' + obj.tangent_n{2}*obj.e2_n')';
322 323
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')'; 324 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')'; 325 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')'; 326 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')'; 327 obj.tau_n_n = (obj.n_n{1}*obj.tau1_n' + obj.n_n{2}*obj.tau2_n')';
328
329 obj.tau_t_w = (obj.tangent_w{1}*obj.tau1_w' + obj.tangent_w{2}*obj.tau2_w')';
330 obj.tau_t_e = (obj.tangent_e{1}*obj.tau1_e' + obj.tangent_e{2}*obj.tau2_e')';
331 obj.tau_t_s = (obj.tangent_s{1}*obj.tau1_s' + obj.tangent_s{2}*obj.tau2_s')';
332 obj.tau_t_n = (obj.tangent_n{1}*obj.tau1_n' + obj.tangent_n{2}*obj.tau2_n')';
332 333
333 % Stress operators 334 % Stress operators
334 sigma = cell(dim, dim); 335 sigma = cell(dim, dim);
335 D1 = {obj.Dx, obj.Dy}; 336 D1 = {obj.Dx, obj.Dy};
336 E = obj.E; 337 E = obj.E;
396 case 't' 397 case 't'
397 en = obj.getBoundaryOperator('et', boundary); 398 en = obj.getBoundaryOperator('et', boundary);
398 end 399 end
399 e1 = obj.getBoundaryOperator('e1', boundary); 400 e1 = obj.getBoundaryOperator('e1', boundary);
400 401
401 bc = {1, type}; 402 bc1 = {1, type};
402 [c1, p1] = obj.refObj.boundary_condition(boundary, bc, tuning); 403 [c1, p1] = obj.refObj.boundary_condition(boundary, bc1, tuning);
403 bc = {2, type}; 404 bc2 = {2, type};
404 c2 = obj.refObj.boundary_condition(boundary, bc, tuning); 405 c2 = obj.refObj.boundary_condition(boundary, bc2, tuning);
405 406
406 switch type 407 switch type
407 case {'F','f','Free','free','traction','Traction','t','T'} 408 case {'F','f','Free','free','traction','Traction','t','T'}
408 closure = en*en'*(c1+c2); 409 closure = en*en'*(c1+c2);
409 penalty = en*e1'*p1; 410 penalty = en*e1'*p1;
410 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'} 411 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
411 error('Not implemented'); 412 [closure, penalty] = obj.displacementBCNormalTangential(boundary, bc, tuning);
412 end 413 end
413 414
414 end 415 end
415 416
416 switch type 417 switch type
421 penalty = penalty*s; 422 penalty = penalty*s;
422 423
423 end 424 end
424 end 425 end
425 426
426 function [closure, penalty] = displacementBCNormalTangential(boundary, bc, tuning) 427 function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning)
427 error('not implemented'); 428 disp('WARNING: displacementBCNormalTangential is only guaranteed to work for displaacement BC on one component and traction on the other');
428 u = obj; 429 u = obj;
429 430
430 component = bc{1}; 431 component = bc{1};
431 type = bc{2}; 432 type = bc{2};
432 433
433 switch component 434 switch component
434 case 'n' 435 case 'n'
435 en_u = u.getBoundaryOperator('en', boundary); 436 en = u.getBoundaryOperator('en', boundary);
436 tau_n_u = u.getBoundaryOperator('tau_n', boundary); 437 tau_n = u.getBoundaryOperator('tau_n', boundary);
438 N = u.getNormal(boundary);
437 case 't' 439 case 't'
438 en_u = u.getBoundaryOperator('et', boundary); 440 en = u.getBoundaryOperator('et', boundary);
439 tau_n_u = u.getBoundaryOperator('tau_t', boundary); 441 tau_n = u.getBoundaryOperator('tau_t', boundary);
442 N = u.getTangent(boundary);
440 end 443 end
441 444
442 % Operators 445 % Operators
446 e = u.getBoundaryOperatorForScalarField('e', boundary);
447 h11 = u.getBorrowing(boundary);
448 n = u.getNormal(boundary);
449
450 C = u.C;
451 Ji = u.Ji;
452 s = spdiag(u.(['s_' boundary]));
453 m_tot = u.grid.N();
454
455 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
456 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
457
458 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
459 dim = u.dim;
460
461 % Preallocate
462 [~, m_int] = size(H_gamma);
463 closure = sparse(dim*m_tot, dim*m_tot);
464 penalty = sparse(dim*m_tot, m_int);
465
466 % Term 1: The symmetric term
467 Z = sparse(m_int, m_int);
468 for i = 1:dim
469 for j = 1:dim
470 for l = 1:dim
471 for k = 1:dim
472 Z = Z + n{i}*N{j}*e'*Ji*C{i,j,k,l}*e*n{k}*N{l};
473 end
474 end
475 end
476 end
477
478 Z = -tuning*dim*1/h11*s*Z;
479 closure = closure + en*H_gamma*Z*en';
480 penalty = penalty - en*H_gamma*Z;
481
482 % Term 2: The symmetrizing term
483 closure = closure + tau_n*H_gamma*en';
484 penalty = penalty - tau_n*H_gamma;
485
486 % Multiply all normal component terms by inverse of density x quadraure
487 closure = RHOi*Hi*closure;
488 penalty = RHOi*Hi*penalty;
489 end
490
491 % type Struct that specifies the interface coupling.
492 % Fields:
493 % -- tuning: penalty strength, defaults to 1.0
494 % -- interpolation: type of interpolation, default 'none'
495 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
496
497 defaultType.tuning = 1.0;
498 defaultType.interpolation = 'none';
499 defaultType.type = 'standard';
500 default_struct('type', defaultType);
501
502 switch type.type
503 case 'standard'
504 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
505 case 'frictionalFault'
506 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
507 end
508
509 end
510
511 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
512 tuning = type.tuning;
513
514 % u denotes the solution in the own domain
515 % v denotes the solution in the neighbour domain
516
517 u = obj;
518 v = neighbour_scheme;
519
520 % Operators, u side
443 e_u = u.getBoundaryOperatorForScalarField('e', boundary); 521 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
522 en_u = u.getBoundaryOperator('en', boundary);
523 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
444 h11_u = u.getBorrowing(boundary); 524 h11_u = u.getBorrowing(boundary);
445 n_u = u.getNormal(boundary); 525 n_u = u.getNormal(boundary);
446 526
447 C_u = u.C; 527 C_u = u.C;
448 Ji_u = u.Ji; 528 Ji_u = u.Ji;
449 s_u = spdiag(u.(['s_' boundary])); 529 s_u = spdiag(u.(['s_' boundary]));
450 m_tot_u = u.grid.N(); 530 m_tot_u = u.grid.N();
451 531
532 % Operators, v side
533 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
534 en_v = v.getBoundaryOperator('en', neighbour_boundary);
535 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
536 h11_v = v.getBorrowing(neighbour_boundary);
537 n_v = v.getNormal(neighbour_boundary);
538
539 C_v = v.C;
540 Ji_v = v.Ji;
541 s_v = spdiag(v.(['s_' neighbour_boundary]));
542 m_tot_v = v.grid.N();
543
544 % Operators that are only required for own domain
452 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}'; 545 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}'; 546 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
454 547
455 % Shared operators 548 % Shared operators
456 H_gamma = u.getBoundaryQuadratureForScalarField(boundary); 549 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
466 Z_v = sparse(m_int, m_int); 559 Z_v = sparse(m_int, m_int);
467 for i = 1:dim 560 for i = 1:dim
468 for j = 1:dim 561 for j = 1:dim
469 for l = 1:dim 562 for l = 1:dim
470 for k = 1:dim 563 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}; 564 Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*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}; 565 Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*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
495 % type Struct that specifies the interface coupling.
496 % Fields:
497 % -- tuning: penalty strength, defaults to 1.0
498 % -- interpolation: type of interpolation, default 'none'
499 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
500
501 defaultType.tuning = 1.0;
502 defaultType.interpolation = 'none';
503 defaultType.type = 'standard';
504 default_struct('type', defaultType);
505
506 switch type.type
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 566 end
571 end 567 end
572 end 568 end
573 end 569 end
574 570
613 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 609 assertIsMember(boundary, {'w', 'e', 's', 'n'})
614 610
615 n = obj.(['n_' boundary]); 611 n = obj.(['n_' boundary]);
616 end 612 end
617 613
614 % Returns the unit tangent vector for the boundary specified by the string boundary.
615 % t is a cell of diagonal matrices for each normal component, t{1} = t_1, t{2} = t_2.
616 function t = getTangent(obj, boundary)
617 assertIsMember(boundary, {'w', 'e', 's', 'n'})
618
619 t = obj.(['tangent_' boundary]);
620 end
621
618 % Returns the boundary operator op for the boundary specified by the string boundary. 622 % Returns the boundary operator op for the boundary specified by the string boundary.
619 % op -- string 623 % op -- string
620 function o = getBoundaryOperator(obj, op, boundary) 624 function o = getBoundaryOperator(obj, op, boundary)
621 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 625 assertIsMember(boundary, {'w', 'e', 's', 'n'})
622 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n'}) 626 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
623 627
624 o = obj.([op, '_', boundary]); 628 o = obj.([op, '_', boundary]);
625 629
626 end 630 end
627 631