Mercurial > repos > public > sbplib
comparison +sbp/+implementations/d2_variable_2.m @ 860:b758d1cf4c8e feature/poroelastic
Add computation of HI*M to D2Variable to make adjoint gradient computation easier.
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Wed, 24 Oct 2018 16:16:43 -0700 |
parents | ded4156e53e2 |
children | e54c2f54dbfe |
comparison
equal
deleted
inserted
replaced
859:4c7532db42cd | 860:b758d1cf4c8e |
---|---|
25 | 25 |
26 % D1 operator | 26 % D1 operator |
27 diags = -1:1; | 27 diags = -1:1; |
28 stencil = [-1/2 0 1/2]; | 28 stencil = [-1/2 0 1/2]; |
29 D1 = stripeMatrix(stencil, diags, m); | 29 D1 = stripeMatrix(stencil, diags, m); |
30 | 30 |
31 D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1; | 31 D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1; |
32 D1(m,m-1)=-1;D1(m,m)=1; | 32 D1(m,m-1)=-1;D1(m,m)=1; |
33 D1=D1/h; | 33 D1=D1/h; |
34 %Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m'); | 34 %Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m'); |
35 | 35 |
38 | 38 |
39 scheme_width = 3; | 39 scheme_width = 3; |
40 scheme_radius = (scheme_width-1)/2; | 40 scheme_radius = (scheme_width-1)/2; |
41 r = (1+scheme_radius):(m-scheme_radius); | 41 r = (1+scheme_radius):(m-scheme_radius); |
42 | 42 |
43 function D2 = D2_fun(c) | 43 function [D2, B] = D2_fun(c) |
44 | 44 |
45 Mm1 = -c(r-1)/2 - c(r)/2; | 45 Mm1 = -c(r-1)/2 - c(r)/2; |
46 M0 = c(r-1)/2 + c(r) + c(r+1)/2; | 46 M0 = c(r-1)/2 + c(r) + c(r+1)/2; |
47 Mp1 = -c(r)/2 - c(r+1)/2; | 47 Mp1 = -c(r)/2 - c(r+1)/2; |
48 | 48 |
52 M(1:2,1:2)=[c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;]; | 52 M(1:2,1:2)=[c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;]; |
53 M(m-1:m,m-1:m)=[c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;]; | 53 M(m-1:m,m-1:m)=[c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;]; |
54 M=M/h; | 54 M=M/h; |
55 | 55 |
56 D2=HI*(-M-c(1)*e_l*d1_l'+c(m)*e_r*d1_r'); | 56 D2=HI*(-M-c(1)*e_l*d1_l'+c(m)*e_r*d1_r'); |
57 B = HI*M; | |
57 end | 58 end |
58 D2 = @D2_fun; | 59 D2 = @D2_fun; |
60 | |
59 end | 61 end |