Mercurial > repos > public > sbplib
annotate +scheme/Elastic2dVariable.m @ 961:2efeedf8c34f feature/poroelastic
Make ElasticVariable2d.boundary_condition use get_boundary_operator for penalty strength.
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Mon, 17 Dec 2018 21:07:59 -0800 |
parents | ac566f3dc9b3 |
children | c75ddd568fcc |
rev | line source |
---|---|
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
1 classdef Elastic2dVariable < scheme.Scheme |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
2 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
3 % Discretizes the elastic wave equation: |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
4 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
5 % opSet should be cell array of opSets, one per dimension. This |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
6 % is useful if we have periodic BC in one direction. |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
7 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
8 properties |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
9 m % Number of points in each direction, possibly a vector |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
10 h % Grid spacing |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
11 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
12 grid |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
13 dim |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
14 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
15 order % Order of accuracy for the approximation |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
16 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
17 % Diagonal matrices for varible coefficients |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
18 LAMBDA % Variable coefficient, related to dilation |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
19 MU % Shear modulus, variable coefficient |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
20 RHO, RHOi % Density, variable |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
21 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
22 D % Total operator |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
23 D1 % First derivatives |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
24 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
25 % Second derivatives |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
26 D2_lambda |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
27 D2_mu |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
28 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
29 % Traction operators used for BC |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
30 T_l, T_r |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
31 tau_l, tau_r |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
32 |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
33 H, Hi, H_1D % Inner products |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
34 e_l, e_r |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
35 d1_l, d1_r % Normal derivatives at the boundary |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
36 E % E{i}^T picks out component i |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
37 |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
38 H_boundary % Boundary inner products |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
39 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
40 % Kroneckered norms and coefficients |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
41 RHOi_kron |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
42 Hi_kron |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
43 |
882
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
44 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant. |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
45 theta_R % Borrowing (d1- D1)^2 from R |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
46 theta_H % First entry in norm matrix |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
47 theta_M % Borrowing d1^2 from M. |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
48 |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
49 % Structures used for adjoint optimization |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
50 B |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
51 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
52 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
53 methods |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
54 |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
55 % The coefficients can either be function handles or grid functions |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
56 function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet) |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
57 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
58 default_arg('lambda', @(x,y) 0*x+1); |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
59 default_arg('mu', @(x,y) 0*x+1); |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
60 default_arg('rho', @(x,y) 0*x+1); |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
61 dim = 2; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
62 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
63 assert(isa(g, 'grid.Cartesian')) |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
64 |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
65 if isa(lambda, 'function_handle') |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
66 lambda = grid.evalOn(g, lambda); |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
67 end |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
68 if isa(mu, 'function_handle') |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
69 mu = grid.evalOn(g, mu); |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
70 end |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
71 if isa(rho, 'function_handle') |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
72 rho = grid.evalOn(g, rho); |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
73 end |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
74 |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
75 m = g.size(); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
76 m_tot = g.N(); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
77 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
78 h = g.scaling(); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
79 lim = g.lim; |
729
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
80 if isempty(lim) |
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
81 x = g.x; |
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
82 lim = cell(length(x),1); |
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
83 for i = 1:length(x) |
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
84 lim{i} = {min(x{i}), max(x{i})}; |
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
85 end |
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
86 end |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
87 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
88 % 1D operators |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
89 ops = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
90 for i = 1:dim |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
91 ops{i} = opSet{i}(m(i), lim{i}, order); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
92 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
93 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
94 % Borrowing constants |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
95 for i = 1:dim |
882
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
96 obj.theta_R{i} = h(i)*ops{i}.borrowing.R.delta_D; |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
97 obj.theta_H{i} = h(i)*ops{i}.borrowing.H11; |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
98 obj.theta_M{i} = h(i)*ops{i}.borrowing.M.d1; |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
99 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
100 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
101 I = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
102 D1 = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
103 D2 = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
104 H = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
105 Hi = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
106 e_l = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
107 e_r = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
108 d1_l = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
109 d1_r = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
110 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
111 for i = 1:dim |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
112 I{i} = speye(m(i)); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
113 D1{i} = ops{i}.D1; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
114 D2{i} = ops{i}.D2; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
115 H{i} = ops{i}.H; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
116 Hi{i} = ops{i}.HI; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
117 e_l{i} = ops{i}.e_l; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
118 e_r{i} = ops{i}.e_r; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
119 d1_l{i} = ops{i}.d1_l; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
120 d1_r{i} = ops{i}.d1_r; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
121 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
122 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
123 %====== Assemble full operators ======== |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
124 LAMBDA = spdiag(lambda); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
125 obj.LAMBDA = LAMBDA; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
126 MU = spdiag(mu); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
127 obj.MU = MU; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
128 RHO = spdiag(rho); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
129 obj.RHO = RHO; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
130 obj.RHOi = inv(RHO); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
131 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
132 obj.D1 = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
133 obj.D2_lambda = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
134 obj.D2_mu = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
135 obj.e_l = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
136 obj.e_r = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
137 obj.d1_l = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
138 obj.d1_r = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
139 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
140 % D1 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
141 obj.D1{1} = kron(D1{1},I{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
142 obj.D1{2} = kron(I{1},D1{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
143 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
144 % Boundary operators |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
145 obj.e_l{1} = kron(e_l{1},I{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
146 obj.e_l{2} = kron(I{1},e_l{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
147 obj.e_r{1} = kron(e_r{1},I{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
148 obj.e_r{2} = kron(I{1},e_r{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
149 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
150 obj.d1_l{1} = kron(d1_l{1},I{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
151 obj.d1_l{2} = kron(I{1},d1_l{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
152 obj.d1_r{1} = kron(d1_r{1},I{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
153 obj.d1_r{2} = kron(I{1},d1_r{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
154 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
155 % D2 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
156 for i = 1:dim |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
157 obj.D2_lambda{i} = sparse(m_tot); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
158 obj.D2_mu{i} = sparse(m_tot); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
159 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
160 ind = grid.funcToMatrix(g, 1:m_tot); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
161 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
162 for i = 1:m(2) |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
163 D_lambda = D2{1}(lambda(ind(:,i))); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
164 D_mu = D2{1}(mu(ind(:,i))); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
165 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
166 p = ind(:,i); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
167 obj.D2_lambda{1}(p,p) = D_lambda; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
168 obj.D2_mu{1}(p,p) = D_mu; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
169 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
170 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
171 for i = 1:m(1) |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
172 D_lambda = D2{2}(lambda(ind(i,:))); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
173 D_mu = D2{2}(mu(ind(i,:))); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
174 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
175 p = ind(i,:); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
176 obj.D2_lambda{2}(p,p) = D_lambda; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
177 obj.D2_mu{2}(p,p) = D_mu; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
178 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
179 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
180 % Quadratures |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
181 obj.H = kron(H{1},H{2}); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
182 obj.Hi = inv(obj.H); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
183 obj.H_boundary = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
184 obj.H_boundary{1} = H{2}; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
185 obj.H_boundary{2} = H{1}; |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
186 obj.H_1D = {H{1}, H{2}}; |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
187 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
188 % E{i}^T picks out component i. |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
189 E = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
190 I = speye(m_tot,m_tot); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
191 for i = 1:dim |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
192 e = sparse(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
193 e(i) = 1; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
194 E{i} = kron(I,e); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
195 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
196 obj.E = E; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
197 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
198 % Differentiation matrix D (without SAT) |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
199 D2_lambda = obj.D2_lambda; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
200 D2_mu = obj.D2_mu; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
201 D1 = obj.D1; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
202 D = sparse(dim*m_tot,dim*m_tot); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
203 d = @kroneckerDelta; % Kronecker delta |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
204 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
205 for i = 1:dim |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
206 for j = 1:dim |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
207 D = D + E{i}*inv(RHO)*( d(i,j)*D2_lambda{i}*E{j}' +... |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
208 db(i,j)*D1{i}*LAMBDA*D1{j}*E{j}' ... |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
209 ); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
210 D = D + E{i}*inv(RHO)*( d(i,j)*D2_mu{i}*E{j}' +... |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
211 db(i,j)*D1{j}*MU*D1{i}*E{j}' + ... |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
212 D2_mu{j}*E{i}' ... |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
213 ); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
214 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
215 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
216 obj.D = D; |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
217 %=========================================%' |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
218 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
219 % Numerical traction operators for BC. |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
220 % Because d1 =/= e0^T*D1, the numerical tractions are different |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
221 % at every boundary. |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
222 T_l = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
223 T_r = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
224 tau_l = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
225 tau_r = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
226 % tau^{j}_i = sum_k T^{j}_{ik} u_k |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
227 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
228 d1_l = obj.d1_l; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
229 d1_r = obj.d1_r; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
230 e_l = obj.e_l; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
231 e_r = obj.e_r; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
232 D1 = obj.D1; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
233 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
234 % Loop over boundaries |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
235 for j = 1:dim |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
236 T_l{j} = cell(dim,dim); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
237 T_r{j} = cell(dim,dim); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
238 tau_l{j} = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
239 tau_r{j} = cell(dim,1); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
240 |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
241 LAMBDA_l = e_l{j}'*LAMBDA*e_l{j}; |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
242 LAMBDA_r = e_r{j}'*LAMBDA*e_r{j}; |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
243 MU_l = e_l{j}'*MU*e_l{j}; |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
244 MU_r = e_r{j}'*MU*e_r{j}; |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
245 |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
246 [~, n_l] = size(e_l{j}); |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
247 [~, n_r] = size(e_r{j}); |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
248 |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
249 % Loop over components |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
250 for i = 1:dim |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
251 tau_l{j}{i} = sparse(n_l, dim*m_tot); |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
252 tau_r{j}{i} = sparse(n_r, dim*m_tot); |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
253 for k = 1:dim |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
254 T_l{j}{i,k} = ... |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
255 -d(i,j)*LAMBDA_l*(d(i,k)*d1_l{k}' + db(i,k)*e_l{j}'*D1{k})... |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
256 -d(j,k)*MU_l*(d(i,j)*d1_l{i}' + db(i,j)*e_l{j}'*D1{i})... |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
257 -d(i,k)*MU_l*d1_l{j}'; |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
258 |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
259 T_r{j}{i,k} = ... |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
260 d(i,j)*LAMBDA_r*(d(i,k)*d1_r{k}' + db(i,k)*e_r{j}'*D1{k})... |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
261 +d(j,k)*MU_r*(d(i,j)*d1_r{i}' + db(i,j)*e_r{j}'*D1{i})... |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
262 +d(i,k)*MU_r*d1_r{j}'; |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
263 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
264 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
265 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
266 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
267 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
268 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
269 end |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
270 |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
271 % Transpose T and tau to match boundary operator convention |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
272 for i = 1:dim |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
273 for j = 1:dim |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
274 tau_l{i}{j} = transpose(tau_l{i}{j}); |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
275 tau_r{i}{j} = transpose(tau_r{i}{j}); |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
276 for k = 1:dim |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
277 T_l{i}{j,k} = transpose(T_l{i}{j,k}); |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
278 T_r{i}{j,k} = transpose(T_r{i}{j,k}); |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
279 end |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
280 end |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
281 end |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
282 |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
283 obj.T_l = T_l; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
284 obj.T_r = T_r; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
285 obj.tau_l = tau_l; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
286 obj.tau_r = tau_r; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
287 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
288 % Kroneckered norms and coefficients |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
289 I_dim = speye(dim); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
290 obj.RHOi_kron = kron(obj.RHOi, I_dim); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
291 obj.Hi_kron = kron(obj.Hi, I_dim); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
292 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
293 % Misc. |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
294 obj.m = m; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
295 obj.h = h; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
296 obj.order = order; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
297 obj.grid = g; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
298 obj.dim = dim; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
299 |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
300 % Used for adjoint optimization |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
301 obj.B = cell(1,dim); |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
302 for i = 1:dim |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
303 obj.B{i} = zeros(m(i),m(i),m(i)); |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
304 for k = 1:m(i) |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
305 c = sparse(m(i),1); |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
306 c(k) = 1; |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
307 [~, obj.B{i}(:,:,k)] = ops{i}.D2(c); |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
308 end |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
309 end |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
310 |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
311 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
312 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
313 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
314 % Closure functions return the operators applied to the own domain to close the boundary |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
315 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
316 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. |
795
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
317 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
318 % on the first component. |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
319 % data is a function returning the data that should be applied at the boundary. |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
320 % neighbour_scheme is an instance of Scheme that should be interfaced to. |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
321 % neighbour_boundary is a string specifying which boundary to interface to. |
795
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
322 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning) |
738
aa4ef495f1fd
Set Dirichlet tuning as parameter in BC function in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
734
diff
changeset
|
323 default_arg('tuning', 1.2); |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
324 |
795
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
325 assert( iscell(bc), 'The BC type must be a 2x1 cell array' ); |
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
326 comp = bc{1}; |
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
327 type = bc{2}; |
729
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
328 |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
329 % j is the coordinate direction of the boundary |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
330 j = obj.get_boundary_number(boundary); |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
331 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
332 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
333 E = obj.E; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
334 Hi = obj.Hi; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
335 LAMBDA = obj.LAMBDA; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
336 MU = obj.MU; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
337 RHOi = obj.RHOi; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
338 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
339 dim = obj.dim; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
340 m_tot = obj.grid.N(); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
341 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
342 % Preallocate |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
343 closure = sparse(dim*m_tot, dim*m_tot); |
795
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
344 penalty = sparse(dim*m_tot, m_tot/obj.m(j)); |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
345 |
795
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
346 k = comp; |
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
347 switch type |
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
348 |
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
349 % Dirichlet boundary condition |
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
350 case {'D','d','dirichlet','Dirichlet'} |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
351 |
961
2efeedf8c34f
Make ElasticVariable2d.boundary_condition use get_boundary_operator for penalty strength.
Martin Almquist <malmquist@stanford.edu>
parents:
960
diff
changeset
|
352 alpha = obj.get_boundary_operator('alpha', boundary); |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
353 |
795
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
354 % Loop over components that Dirichlet penalties end up on |
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
355 for i = 1:dim |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
356 C = transpose(T{k,i}); |
961
2efeedf8c34f
Make ElasticVariable2d.boundary_condition use get_boundary_operator for penalty strength.
Martin Almquist <malmquist@stanford.edu>
parents:
960
diff
changeset
|
357 A = -alpha{i,k}; |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
358 B = A + e*C; |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
359 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); |
795
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
360 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
361 end |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
362 |
795
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
363 % Free boundary condition |
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
364 case {'F','f','Free','free','traction','Traction','t','T'} |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
365 closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau{k}'; |
795
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
366 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
367 |
795
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
368 % Unknown boundary condition |
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
369 otherwise |
1f6b2fb69225
Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents:
738
diff
changeset
|
370 error('No such boundary condition: type = %s',type); |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
371 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
372 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
373 |
960
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
374 % type Struct that specifies the interface coupling. |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
375 % Fields: |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
376 % -- tuning: penalty strength, defaults to 1.2 |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
377 % -- interpolation: type of interpolation, default 'none' |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
378 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
379 |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
380 defaultType.tuning = 1.2; |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
381 defaultType.interpolation = 'none'; |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
382 default_struct('type', defaultType); |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
383 |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
384 switch type.interpolation |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
385 case {'none', ''} |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
386 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
387 case {'op','OP'} |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
388 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
389 otherwise |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
390 error('Unknown type of interpolation: %s ', type.interpolation); |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
391 end |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
392 end |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
393 |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
394 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
395 % u denotes the solution in the own domain |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
396 % v denotes the solution in the neighbour domain |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
397 % Operators without subscripts are from the own domain. |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
398 |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
399 % j is the coordinate direction of the boundary |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
400 j = obj.get_boundary_number(boundary); |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
401 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary); |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
402 |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
403 % Get boundary operators |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
404 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
405 [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary); |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
406 |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
407 % Operators and quantities that correspond to the own domain only |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
408 Hi = obj.Hi; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
409 RHOi = obj.RHOi; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
410 dim = obj.dim; |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
411 |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
412 %--- Other operators ---- |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
413 m_tot_u = obj.grid.N(); |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
414 E = obj.E; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
415 LAMBDA_u = obj.LAMBDA; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
416 MU_u = obj.MU; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
417 lambda_u = e'*LAMBDA_u*e; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
418 mu_u = e'*MU_u*e; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
419 |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
420 m_tot_v = neighbour_scheme.grid.N(); |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
421 E_v = neighbour_scheme.E; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
422 LAMBDA_v = neighbour_scheme.LAMBDA; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
423 MU_v = neighbour_scheme.MU; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
424 lambda_v = e_v'*LAMBDA_v*e_v; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
425 mu_v = e_v'*MU_v*e_v; |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
426 %------------------------- |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
427 |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
428 % Borrowing constants |
882
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
429 theta_R_u = obj.theta_R{j}; |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
430 theta_H_u = obj.theta_H{j}; |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
431 theta_M_u = obj.theta_M{j}; |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
432 |
882
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
433 theta_R_v = neighbour_scheme.theta_R{j_v}; |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
434 theta_H_v = neighbour_scheme.theta_H{j_v}; |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
435 theta_M_v = neighbour_scheme.theta_M{j_v}; |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
436 |
882
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
437 function [alpha_ii, alpha_ij] = computeAlpha(th_R, th_H, th_M, lambda, mu) |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
438 alpha_ii = dim*lambda/(4*th_H) + lambda/(4*th_R) + mu/(2*th_M); |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
439 alpha_ij = mu/(2*th_H) + mu/(4*th_R); |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
440 end |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
441 |
882
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
442 [alpha_ii_u, alpha_ij_u] = computeAlpha(theta_R_u, theta_H_u, theta_M_u, lambda_u, mu_u); |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
443 [alpha_ii_v, alpha_ij_v] = computeAlpha(theta_R_v, theta_H_v, theta_M_v, lambda_v, mu_v); |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
444 sigma_ii = tuning*(alpha_ii_u + alpha_ii_v); |
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
445 sigma_ij = tuning*(alpha_ij_u + alpha_ij_v); |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
446 |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
447 d = @kroneckerDelta; % Kronecker delta |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
448 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
449 sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij); |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
450 |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
451 % Preallocate |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
452 closure = sparse(dim*m_tot_u, dim*m_tot_u); |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
453 penalty = sparse(dim*m_tot_u, dim*m_tot_v); |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
454 |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
455 % Loop over components that penalties end up on |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
456 for i = 1:dim |
729
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
457 closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}'; |
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
458 penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}'; |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
459 |
729
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
460 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i}; |
aa8cf3851de8
Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents:
727
diff
changeset
|
461 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i}; |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
462 |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
463 % Loop over components that we have interface conditions on |
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
464 for k = 1:dim |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
465 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e'*E{k}'; |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
466 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e_v'*E_v{k}'; |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
467 end |
727
6d5953fc090e
First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents:
726
diff
changeset
|
468 end |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
469 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
470 |
960
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
471 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
472 error('Non-conforming interfaces not implemented yet.'); |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
473 end |
ac566f3dc9b3
Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents:
959
diff
changeset
|
474 |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
475 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
476 function [j, nj] = get_boundary_number(obj, boundary) |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
477 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
478 switch boundary |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
479 case {'w','W','west','West', 'e', 'E', 'east', 'East'} |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
480 j = 1; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
481 case {'s','S','south','South', 'n', 'N', 'north', 'North'} |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
482 j = 2; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
483 otherwise |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
484 error('No such boundary: boundary = %s',boundary); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
485 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
486 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
487 switch boundary |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
488 case {'w','W','west','West','s','S','south','South'} |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
489 nj = -1; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
490 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
491 nj = 1; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
492 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
493 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
494 |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
495 % Returns the boundary operator op for the boundary specified by the string boundary. |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
496 % op: may be a cell array of strings |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
497 function [varargout] = get_boundary_operator(obj, op, boundary) |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
498 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
499 switch boundary |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
500 case {'w','W','west','West', 'e', 'E', 'east', 'East'} |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
501 j = 1; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
502 case {'s','S','south','South', 'n', 'N', 'north', 'North'} |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
503 j = 2; |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
504 otherwise |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
505 error('No such boundary: boundary = %s',boundary); |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
506 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
507 |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
508 if ~iscell(op) |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
509 op = {op}; |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
510 end |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
511 |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
512 for k = 1:length(op) |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
513 switch op{k} |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
514 case 'e' |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
515 switch boundary |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
516 case {'w','W','west','West','s','S','south','South'} |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
517 varargout{k} = obj.e_l{j}; |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
518 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
519 varargout{k} = obj.e_r{j}; |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
520 end |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
521 case 'd' |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
522 switch boundary |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
523 case {'w','W','west','West','s','S','south','South'} |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
524 varargout{k} = obj.d1_l{j}; |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
525 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
526 varargout{k} = obj.d1_r{j}; |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
527 end |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
528 case 'H' |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
529 varargout{k} = obj.H_boundary{j}; |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
530 case 'T' |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
531 switch boundary |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
532 case {'w','W','west','West','s','S','south','South'} |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
533 varargout{k} = obj.T_l{j}; |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
534 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
535 varargout{k} = obj.T_r{j}; |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
536 end |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
537 case 'tau' |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
538 switch boundary |
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
539 case {'w','W','west','West','s','S','south','South'} |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
540 varargout{k} = obj.tau_l{j}; |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
541 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
542 varargout{k} = obj.tau_r{j}; |
855
5751262b323b
Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents:
795
diff
changeset
|
543 end |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
544 case 'alpha' |
882
14fee299ada2
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
Martin Almquist <malmquist@stanford.edu>
parents:
861
diff
changeset
|
545 % alpha = alpha(i,j) is the penalty strength for displacement BC. |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
546 tuning = 1.2; |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
547 LAMBDA = obj.LAMBDA; |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
548 MU = obj.MU; |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
549 |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
550 dim = obj.dim; |
919
e30aaa4a3e09
Bugfix in Elastic2dVariable.get_boundary_ops
Martin Almquist <malmquist@stanford.edu>
parents:
882
diff
changeset
|
551 theta_R = obj.theta_R{j}; |
e30aaa4a3e09
Bugfix in Elastic2dVariable.get_boundary_ops
Martin Almquist <malmquist@stanford.edu>
parents:
882
diff
changeset
|
552 theta_H = obj.theta_H{j}; |
e30aaa4a3e09
Bugfix in Elastic2dVariable.get_boundary_ops
Martin Almquist <malmquist@stanford.edu>
parents:
882
diff
changeset
|
553 theta_M = obj.theta_M{j}; |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
554 |
919
e30aaa4a3e09
Bugfix in Elastic2dVariable.get_boundary_ops
Martin Almquist <malmquist@stanford.edu>
parents:
882
diff
changeset
|
555 a_lambda = dim/theta_H + 1/theta_R; |
e30aaa4a3e09
Bugfix in Elastic2dVariable.get_boundary_ops
Martin Almquist <malmquist@stanford.edu>
parents:
882
diff
changeset
|
556 a_mu_i = 2/theta_M; |
e30aaa4a3e09
Bugfix in Elastic2dVariable.get_boundary_ops
Martin Almquist <malmquist@stanford.edu>
parents:
882
diff
changeset
|
557 a_mu_ij = 2/theta_H + 1/theta_R; |
861
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
558 |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
559 d = @kroneckerDelta; % Kronecker delta |
607c631f175e
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents:
855
diff
changeset
|
560 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
561 alpha = cell(obj.dim, obj.dim); |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
562 |
959
c226fb8c2b8a
Bugfix in Elastic2dVariable.getBoundaryOps, alpha. Adjoint FD conv works now!
Martin Almquist <malmquist@stanford.edu>
parents:
958
diff
changeset
|
563 alpha_func = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... |
919
e30aaa4a3e09
Bugfix in Elastic2dVariable.get_boundary_ops
Martin Almquist <malmquist@stanford.edu>
parents:
882
diff
changeset
|
564 + d(i,j)* a_mu_i*MU ... |
e30aaa4a3e09
Bugfix in Elastic2dVariable.get_boundary_ops
Martin Almquist <malmquist@stanford.edu>
parents:
882
diff
changeset
|
565 + db(i,j)*a_mu_ij*MU ); |
959
c226fb8c2b8a
Bugfix in Elastic2dVariable.getBoundaryOps, alpha. Adjoint FD conv works now!
Martin Almquist <malmquist@stanford.edu>
parents:
958
diff
changeset
|
566 for i = 1:obj.dim |
c226fb8c2b8a
Bugfix in Elastic2dVariable.getBoundaryOps, alpha. Adjoint FD conv works now!
Martin Almquist <malmquist@stanford.edu>
parents:
958
diff
changeset
|
567 for l = 1:obj.dim |
c226fb8c2b8a
Bugfix in Elastic2dVariable.getBoundaryOps, alpha. Adjoint FD conv works now!
Martin Almquist <malmquist@stanford.edu>
parents:
958
diff
changeset
|
568 alpha{i,l} = d(i,l)*alpha_func(i,j); |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
569 end |
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
570 end |
919
e30aaa4a3e09
Bugfix in Elastic2dVariable.get_boundary_ops
Martin Almquist <malmquist@stanford.edu>
parents:
882
diff
changeset
|
571 |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
572 varargout{k} = alpha; |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
573 otherwise |
958
72cd29107a9a
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
Martin Almquist <malmquist@stanford.edu>
parents:
919
diff
changeset
|
574 error(['No such operator: operator = ' op{k}]); |
726
37d5d69b1a0d
Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents:
689
diff
changeset
|
575 end |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
576 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
577 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
578 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
579 |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
580 function N = size(obj) |
734
eebe24a636c7
Make Elastic2dVariable.size account for components.
Martin Almquist <malmquist@stanford.edu>
parents:
730
diff
changeset
|
581 N = obj.dim*prod(obj.m); |
687
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
582 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
583 end |
e8fc3aa1faf6
Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
584 end |