Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dStaggeredAnisotropic.m @ 1265:6696b216b982 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Start working on displacement bc, clean up traction.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 30 Apr 2020 11:30:55 -0700
parents 066fdfaa2411
children ad31d9c4cec2
comparison
equal deleted inserted replaced
1264:066fdfaa2411 1265:6696b216b982
28 % Boundary operators in cell format, used for BC 28 % Boundary operators in cell format, used for BC
29 % T_w, T_e, T_s, T_n 29 % T_w, T_e, T_s, T_n
30 30
31 % Traction operators 31 % Traction operators
32 tau_w, tau_e, tau_s, tau_n % Return scalar field 32 tau_w, tau_e, tau_s, tau_n % Return scalar field
33 T_w, T_e, T_s, T_n % Act on scalar, return scalar
33 34
34 % Inner products 35 % Inner products
35 H, H_u, H_s 36 H, H_u, H_s
36 % , Hi, Hi_kron, H_1D 37 % , Hi, Hi_kron, H_1D
37 38
343 tau_w = cell(nGrids, 1); 344 tau_w = cell(nGrids, 1);
344 tau_e = cell(nGrids, 1); 345 tau_e = cell(nGrids, 1);
345 tau_s = cell(nGrids, 1); 346 tau_s = cell(nGrids, 1);
346 tau_n = cell(nGrids, 1); 347 tau_n = cell(nGrids, 1);
347 348
349 T_w = cell(nGrids, nGrids);
350 T_e = cell(nGrids, nGrids);
351 T_s = cell(nGrids, nGrids);
352 T_n = cell(nGrids, nGrids);
353 for b = 1:nGrids
354 [~, m_we] = size(e_w_s{b});
355 [~, m_sn] = size(e_s_s{b});
356 for c = 1:nGrids
357 T_w{b,c} = cell(dim, dim);
358 T_e{b,c} = cell(dim, dim);
359 T_s{b,c} = cell(dim, dim);
360 T_n{b,c} = cell(dim, dim);
361
362 for i = 1:dim
363 for j = 1:dim
364 T_w{b,c}{i,j} = sparse(m_we, N_u{c});
365 T_e{b,c}{i,j} = sparse(m_we, N_u{c});
366 T_s{b,c}{i,j} = sparse(m_sn, N_u{c});
367 T_n{b,c}{i,j} = sparse(m_sn, N_u{c});
368 end
369 end
370 end
371 end
372
348 for b = 1:nGrids 373 for b = 1:nGrids
349 tau_w{b} = cell(dim, 1); 374 tau_w{b} = cell(dim, 1);
350 tau_e{b} = cell(dim, 1); 375 tau_e{b} = cell(dim, 1);
351 tau_s{b} = cell(dim, 1); 376 tau_s{b} = cell(dim, 1);
352 tau_n{b} = cell(dim, 1); 377 tau_n{b} = cell(dim, 1);
353 378
354 for j = 1:dim 379 for j = 1:dim
355 tau_w{b}{j} = sparse(N, m_s{b}(2)); 380 tau_w{b}{j} = sparse(N, m_s{b}(2));
356 tau_e{b}{j} = sparse(N, m_s{b}(2)); 381 tau_e{b}{j} = sparse(N, m_s{b}(2));
357 tau_s{b}{j} = sparse(N, m_s{b}(1)); 382 tau_s{b}{j} = sparse(N, m_s{b}(1));
358 tau_n{b}{j} = sparse(N, m_s{b}(1)); 383 tau_n{b}{j} = sparse(N, m_s{b}(1));
367 392
368 tau_w{b}{j} = tau_w{b}{j} + (e_w_s{b}'*n_w(i)*sigma_b_ij)'; 393 tau_w{b}{j} = tau_w{b}{j} + (e_w_s{b}'*n_w(i)*sigma_b_ij)';
369 tau_e{b}{j} = tau_e{b}{j} + (e_e_s{b}'*n_e(i)*sigma_b_ij)'; 394 tau_e{b}{j} = tau_e{b}{j} + (e_e_s{b}'*n_e(i)*sigma_b_ij)';
370 tau_s{b}{j} = tau_s{b}{j} + (e_s_s{b}'*n_s(i)*sigma_b_ij)'; 395 tau_s{b}{j} = tau_s{b}{j} + (e_s_s{b}'*n_s(i)*sigma_b_ij)';
371 tau_n{b}{j} = tau_n{b}{j} + (e_n_s{b}'*n_n(i)*sigma_b_ij)'; 396 tau_n{b}{j} = tau_n{b}{j} + (e_n_s{b}'*n_n(i)*sigma_b_ij)';
397
398 S_bc_ijl = C{b}{i,j,k,l}*D1_u2s{b,c}{k};
399
400 % T_w{b,c}{i,l} = T_w{b,c}{i,l} + (e_w_s{b}'*n_w(i)*S_bc_ijl)';
401 % T_e{b,c}{i,l} = T_e{b,c}{i,l} + (e_e_s{b}'*n_e(i)*S_bc_ijl)';
402 % T_s{b,c}{i,l} = T_s{b,c}{i,l} + (e_s_s{b}'*n_s(i)*S_bc_ijl)';
403 % T_n{b,c}{i,l} = T_n{b,c}{i,l} + (e_n_s{b}'*n_n(i)*S_bc_ijl)';
372 end 404 end
373 end 405 end
374 end 406 end
375 end 407 end
376 end 408 end
378 410
379 obj.tau_w = tau_w; 411 obj.tau_w = tau_w;
380 obj.tau_e = tau_e; 412 obj.tau_e = tau_e;
381 obj.tau_s = tau_s; 413 obj.tau_s = tau_s;
382 obj.tau_n = tau_n; 414 obj.tau_n = tau_n;
415
416 obj.T_w = T_w;
417 obj.T_e = T_e;
418 obj.T_s = T_s;
419 obj.T_n = T_n;
383 420
384 % D1 = obj.D1; 421 % D1 = obj.D1;
385 422
386 % Traction tensors, T_ij 423 % Traction tensors, T_ij
387 % obj.T_w = T_l{1}; 424 % obj.T_w = T_l{1};
463 type = bc{2}; 500 type = bc{2};
464 if ischar(comp) 501 if ischar(comp)
465 comp = obj.getComponent(comp, boundary); 502 comp = obj.getComponent(comp, boundary);
466 end 503 end
467 504
468 e = obj.getBoundaryOperatorForScalarField('e_u', boundary); 505 e_u = obj.getBoundaryOperatorForScalarField('e_u', boundary);
506 e_s = obj.getBoundaryOperatorForScalarField('e_s', boundary);
469 tau = obj.getBoundaryOperator('tau', boundary); 507 tau = obj.getBoundaryOperator('tau', boundary);
470 % T = obj.getBoundaryTractionOperator(boundary); 508 T = obj.getBoundaryTractionOperator(boundary);
471 % h11 = obj.getBorrowing(boundary); 509 % h11 = obj.getBorrowing(boundary);
472 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary); 510 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
473 nu = obj.getNormal(boundary); 511 nu = obj.getNormal(boundary);
474 512
475 U = obj.U; 513 U = obj.U;
476 G = obj.G; 514 G = obj.G;
477 H = obj.H_u; 515 H = obj.H_u;
478 RHO = obj.RHO; 516 RHO = obj.RHO;
479 C = obj.C; 517 C = obj.C;
480 518
519 %---- Grid layout -------
520 % gu1 = xp o yp;
521 % gu2 = xd o yd;
522 % gs1 = xd o yp;
523 % gs2 = xp o yd;
524 %------------------------
525
481 switch boundary 526 switch boundary
482 case {'s', 'n'} 527 case {'w', 'e'}
483 e = rot90(e,2); 528 gridCombos = {{1,1}, {2,2}};
484 G = rot90(G,2); 529 case {'s', 'n'}
485 U = rot90(U,2); 530 gridCombos = {{2,1}, {1,2}};
486 H = rot90(H,2);
487 RHO = rot90(RHO,2);
488 end 531 end
489 532
490 dim = obj.dim; 533 dim = obj.dim;
491 nGrids = obj.nGrids; 534 nGrids = obj.nGrids;
492 535
493 m_tot = obj.N; 536 m_tot = obj.N;
494 537
495 % Preallocate 538 % Preallocate
496 [~, col] = size(tau{1}{1}); 539 [~, col] = size(tau{1}{1});
497 closure = sparse(m_tot, m_tot); 540 closure = sparse(m_tot, m_tot);
498 % penalty = sparse(m_tot, col);
499 penalty = cell(nGrids, 1); 541 penalty = cell(nGrids, 1);
500 542
501 j = comp; 543 j = comp;
502 switch type 544 switch type
503 545
504 % Dirichlet boundary condition 546 % Dirichlet boundary condition
505 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'} 547 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
506 error('not implemented');
507 548
508 if numel(bc) >= 3 549 if numel(bc) >= 3
509 dComps = bc{3}; 550 dComps = bc{3};
510 else 551 else
511 dComps = 1:dim; 552 dComps = 1:dim;
516 % Z: symmetric part of penalty 557 % Z: symmetric part of penalty
517 % X = Y + Z. 558 % X = Y + Z.
518 559
519 % Nonsymmetric part goes on all components to 560 % Nonsymmetric part goes on all components to
520 % yield traction in discrete energy rate 561 % yield traction in discrete energy rate
521 for i = 1:dim 562 for a = 1:nGrids
522 Y = T{j,i}'; 563 for b = 1:nGrids
523 X = e*Y; 564 for i = 1:dim
524 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' ); 565 Y = T{a,b}{j,i}';
525 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma; 566 X = e_s{a}*Y;
567 closure = closure + U{i}*RHOi*Hi*X'*e_s{a}*H_gamma{a}*(e_u{b}'*U{j}'*G{b}' );
568 penalty = penalty - U{i}*RHOi*Hi*X'*e_s{a}*H_gamma{a};
569 end
570 end
526 end 571 end
527 572
528 % Symmetric part only required on components with displacement BC. 573 % Symmetric part only required on components with displacement BC.
529 % (Otherwise it's not symmetric.) 574 % (Otherwise it's not symmetric.)
530 for i = dComps 575 for i = dComps
540 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma; 585 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
541 end 586 end
542 587
543 % Free boundary condition 588 % Free boundary condition
544 case {'F','f','Free','free','traction','Traction','t','T'} 589 case {'F','f','Free','free','traction','Traction','t','T'}
545 for a = 1:nGrids 590 for m = 1:numel(gridCombos)
546 closure = closure - G{a}*U{a}{j}*(RHO{a}\(H{a}\(e{a}*H_gamma{a}*tau{a}{j}'))); 591 gc = gridCombos{m};
547 penalty{a} = G{a}*U{a}{j}*(RHO{a}\(H{a}\(e{a}*H_gamma{a}))); 592 a = gc{1};
593 b = gc{2};
594 closure = closure - G{a}*U{a}{j}*(RHO{a}\(H{a}\(e_u{a}*H_gamma{b}*tau{b}{j}')));
595 penalty{b} = G{a}*U{a}{j}*(RHO{a}\(H{a}\(e_u{a}*H_gamma{b})));
548 end 596 end
549 597
550 % Unknown boundary condition 598 % Unknown boundary condition
551 otherwise 599 otherwise
552 error('No such boundary condition: type = %s',type); 600 error('No such boundary condition: type = %s',type);