annotate +scheme/Elastic2dVariable.m @ 1208:7f427275bc9c feature/poroelastic

Merge heads
author Martin Almquist <malmquist@stanford.edu>
date Fri, 20 Sep 2019 15:00:19 -0700
parents d9da4c1cdaa0
children 70939ea9a71f
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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:
813
b374a8aa9246 Correct interface penalty strength 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
1200
d9da4c1cdaa0 Fix spelling mistake in comments.
Martin Almquist <malmquist@stanford.edu>
parents: 1118
diff changeset
17 % Diagonal matrices for variable coefficients
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
18 LAMBDA % Lame's first parameter, related to dilation
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
19 MU % Shear modulus
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
20 RHO, RHOi, RHOi_kron % Density
687
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
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
29 % Boundary operators in cell format, used for BC
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
30 T_w, T_e, T_s, T_n
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
32 % Traction operators
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
33 tau_w, tau_e, tau_s, tau_n % Return vector field
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
34 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
35 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
36
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
37 % Inner products
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
38 H, Hi, Hi_kron, H_1D
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
39
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
40 % Boundary inner products (for scalar field)
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
41 H_w, H_e, H_s, H_n
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
42
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
43 % Boundary restriction operators
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
44 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
45 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
46 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
47 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
49 % E{i}^T picks out component i
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
50 E
861
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
51
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
52 % 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
53 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
54 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
55 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
56
861
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
57 % Structures used for adjoint optimization
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
58 B
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 methods
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62
861
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
63 % The coefficients can either be function handles or grid functions
1118
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
64 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
65 function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet, optFlag)
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 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
67 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
68 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
69 default_arg('rho', @(x,y) 0*x+1);
1118
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
70 default_arg('optFlag', false);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 dim = 2;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 assert(isa(g, 'grid.Cartesian'))
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74
861
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
75 if isa(lambda, 'function_handle')
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
76 lambda = grid.evalOn(g, lambda);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
77 end
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
78 if isa(mu, 'function_handle')
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
79 mu = grid.evalOn(g, mu);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
80 end
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
81 if isa(rho, 'function_handle')
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
82 rho = grid.evalOn(g, rho);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
83 end
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
84
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 m = g.size();
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 m_tot = g.N();
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 h = g.scaling();
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 lim = g.lim;
729
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
90 if isempty(lim)
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
91 x = g.x;
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
92 lim = cell(length(x),1);
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
93 for i = 1:length(x)
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
94 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
95 end
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
96 end
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 % 1D operators
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 ops = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101 ops{i} = opSet{i}(m(i), lim{i}, order);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 % Borrowing constants
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 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
106 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
107 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
108 obj.theta_M{i} = h(i)*ops{i}.borrowing.M.d1;
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 end
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 I = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 D1 = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113 D2 = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 H = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 Hi = cell(dim,1);
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
116 e_0 = cell(dim,1);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
117 e_m = cell(dim,1);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
118 d1_0 = cell(dim,1);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
119 d1_m = cell(dim,1);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 I{i} = speye(m(i));
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 D1{i} = ops{i}.D1;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 D2{i} = ops{i}.D2;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 H{i} = ops{i}.H;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 Hi{i} = ops{i}.HI;
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
127 e_0{i} = ops{i}.e_l;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
128 e_m{i} = ops{i}.e_r;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
129 d1_0{i} = ops{i}.d1_l;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
130 d1_m{i} = ops{i}.d1_r;
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 %====== Assemble full operators ========
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134 LAMBDA = spdiag(lambda);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135 obj.LAMBDA = LAMBDA;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 MU = spdiag(mu);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 obj.MU = MU;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 RHO = spdiag(rho);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139 obj.RHO = RHO;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140 obj.RHOi = inv(RHO);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142 obj.D1 = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143 obj.D2_lambda = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144 obj.D2_mu = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 % D1
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 obj.D1{1} = kron(D1{1},I{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148 obj.D1{2} = kron(I{1},D1{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
150 % Boundary restriction operators
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
151 e_l = cell(dim,1);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
152 e_r = cell(dim,1);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
153 e_l{1} = kron(e_0{1}, I{2});
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
154 e_l{2} = kron(I{1}, e_0{2});
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
155 e_r{1} = kron(e_m{1}, I{2});
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
156 e_r{2} = kron(I{1}, e_m{2});
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
157
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
158 e_scalar_w = e_l{1};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
159 e_scalar_e = e_r{1};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
160 e_scalar_s = e_l{2};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
161 e_scalar_n = e_r{2};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
162
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
163 I_dim = speye(dim, dim);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
164 e_w = kron(e_scalar_w, I_dim);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
165 e_e = kron(e_scalar_e, I_dim);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
166 e_s = kron(e_scalar_s, I_dim);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
167 e_n = kron(e_scalar_n, I_dim);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
169 % Boundary derivatives
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
170 d1_l = cell(dim,1);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
171 d1_r = cell(dim,1);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
172 d1_l{1} = kron(d1_0{1}, I{2});
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
173 d1_l{2} = kron(I{1}, d1_0{2});
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
174 d1_r{1} = kron(d1_m{1}, I{2});
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
175 d1_r{2} = kron(I{1}, d1_m{2});
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
176
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
177
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
178 % E{i}^T picks out component i.
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
179 E = cell(dim,1);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
180 I = speye(m_tot,m_tot);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
181 for i = 1:dim
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
182 e = sparse(dim,1);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
183 e(i) = 1;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
184 E{i} = kron(I,e);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
185 end
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
186 obj.E = E;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
187
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
188 e1_w = (e_scalar_w'*E{1}')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
189 e1_e = (e_scalar_e'*E{1}')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
190 e1_s = (e_scalar_s'*E{1}')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
191 e1_n = (e_scalar_n'*E{1}')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
192
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
193 e2_w = (e_scalar_w'*E{2}')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
194 e2_e = (e_scalar_e'*E{2}')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
195 e2_s = (e_scalar_s'*E{2}')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
196 e2_n = (e_scalar_n'*E{2}')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
197
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 % D2
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201 obj.D2_lambda{i} = sparse(m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 obj.D2_mu{i} = sparse(m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 ind = grid.funcToMatrix(g, 1:m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206 for i = 1:m(2)
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207 D_lambda = D2{1}(lambda(ind(:,i)));
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 D_mu = D2{1}(mu(ind(:,i)));
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 p = ind(:,i);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 obj.D2_lambda{1}(p,p) = D_lambda;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212 obj.D2_mu{1}(p,p) = D_mu;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
213 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 for i = 1:m(1)
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 D_lambda = D2{2}(lambda(ind(i,:)));
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 D_mu = D2{2}(mu(ind(i,:)));
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 p = ind(i,:);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 obj.D2_lambda{2}(p,p) = D_lambda;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221 obj.D2_mu{2}(p,p) = D_mu;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 % Quadratures
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 obj.H = kron(H{1},H{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226 obj.Hi = inv(obj.H);
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
227 obj.H_w = H{2};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
228 obj.H_e = H{2};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
229 obj.H_s = H{1};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
230 obj.H_n = H{1};
855
5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable.
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
231 obj.H_1D = {H{1}, H{2}};
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
232
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
233 % Differentiation matrix D (without SAT)
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
234 D2_lambda = obj.D2_lambda;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
235 D2_mu = obj.D2_mu;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
236 D1 = obj.D1;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
237 D = sparse(dim*m_tot,dim*m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
238 d = @kroneckerDelta; % Kronecker delta
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
239 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241 for j = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
242 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
243 db(i,j)*D1{i}*LAMBDA*D1{j}*E{j}' ...
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
244 );
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
245 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
246 db(i,j)*D1{j}*MU*D1{i}*E{j}' + ...
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
247 D2_mu{j}*E{i}' ...
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248 );
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
249 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
250 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
251 obj.D = D;
861
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
252 %=========================================%'
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
254 % Numerical traction operators for BC.
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 % Because d1 =/= e0^T*D1, the numerical tractions are different
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256 % at every boundary.
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
257 %
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
258 % Formula at boundary j: % tau^{j}_i = sum_k T^{j}_{ik} u_k
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
259 %
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
260 T_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
261 T_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
262 tau_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
263 tau_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
264
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
265 D1 = obj.D1;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
266
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
267 % Loop over boundaries
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
268 for j = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
269 T_l{j} = cell(dim,dim);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270 T_r{j} = cell(dim,dim);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
271 tau_l{j} = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
272 tau_r{j} = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
273
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
274 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
275 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
276 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
277 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
278
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 [~, 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
280 [~, 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
281
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
282 % Loop over components
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
283 for i = 1:dim
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
284 tau_l{j}{i} = sparse(dim*m_tot, n_l);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
285 tau_r{j}{i} = sparse(dim*m_tot, n_r);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
286 for k = 1:dim
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
287 T_l{j}{i,k} = ...
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
288 (-d(i,j)*LAMBDA_l*(d(i,k)*d1_l{j}' + db(i,k)*e_l{j}'*D1{k})...
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
289 -d(j,k)*MU_l*(d(i,j)*d1_l{j}' + db(i,j)*e_l{j}'*D1{i})...
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
290 -d(i,k)*MU_l*d1_l{j}')';
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
291
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
292 T_r{j}{i,k} = ...
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
293 (d(i,j)*LAMBDA_r*(d(i,k)*d1_r{j}' + db(i,k)*e_r{j}'*D1{k})...
966
4a85cff53107 Fix bug in Elastic2dVariable that caused crash for m_x neq m_y
Martin Almquist <malmquist@stanford.edu>
parents: 965
diff changeset
294 +d(j,k)*MU_r*(d(i,j)*d1_r{j}' + db(i,j)*e_r{j}'*D1{i})...
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
295 +d(i,k)*MU_r*d1_r{j}')';
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
296
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
297 tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,k}'*E{k}')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
298 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,k}'*E{k}')';
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
299 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
300
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
301 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
302 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
303
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
304 % Traction tensors, T_ij
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
305 obj.T_w = T_l{1};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
306 obj.T_e = T_r{1};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
307 obj.T_s = T_l{2};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
308 obj.T_n = T_r{2};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
309
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
310 % Restriction operators
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
311 obj.e_w = e_w;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
312 obj.e_e = e_e;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
313 obj.e_s = e_s;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
314 obj.e_n = e_n;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
315
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
316 obj.e1_w = e1_w;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
317 obj.e1_e = e1_e;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
318 obj.e1_s = e1_s;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
319 obj.e1_n = e1_n;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
320
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
321 obj.e2_w = e2_w;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
322 obj.e2_e = e2_e;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
323 obj.e2_s = e2_s;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
324 obj.e2_n = e2_n;
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
325
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
326 obj.e_scalar_w = e_scalar_w;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
327 obj.e_scalar_e = e_scalar_e;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
328 obj.e_scalar_s = e_scalar_s;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
329 obj.e_scalar_n = e_scalar_n;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
330
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
331 % First component of traction
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
332 obj.tau1_w = tau_l{1}{1};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
333 obj.tau1_e = tau_r{1}{1};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
334 obj.tau1_s = tau_l{2}{1};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
335 obj.tau1_n = tau_r{2}{1};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
336
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
337 % Second component of traction
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
338 obj.tau2_w = tau_l{1}{2};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
339 obj.tau2_e = tau_r{1}{2};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
340 obj.tau2_s = tau_l{2}{2};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
341 obj.tau2_n = tau_r{2}{2};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
342
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
343 % Traction vectors
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
344 obj.tau_w = (e_w'*e1_w*obj.tau1_w')' + (e_w'*e2_w*obj.tau2_w')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
345 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
346 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')';
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
347 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
348
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
349 % Kroneckered norms and coefficients
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
350 obj.RHOi_kron = kron(obj.RHOi, I_dim);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 obj.Hi_kron = kron(obj.Hi, I_dim);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
353 % Misc.
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
354 obj.m = m;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
355 obj.h = h;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
356 obj.order = order;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
357 obj.grid = g;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
358 obj.dim = dim;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
359
965
db3411264b96 Remove use of tensor library in B assembly
Martin Almquist <malmquist@stanford.edu>
parents: 963
diff changeset
360 % B, used for adjoint optimization
1118
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
361 B = [];
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
362 if optFlag
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
363 B = cell(dim, 1);
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
364 for i = 1:dim
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
365 B{i} = cell(m_tot, 1);
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
366 end
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
367
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
368 B0 = sparse(m_tot, m_tot);
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
369 for i = 1:dim
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
370 for j = 1:m_tot
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
371 B{i}{j} = B0;
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
372 end
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
373 end
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
374
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
375 ind = grid.funcToMatrix(g, 1:m_tot);
965
db3411264b96 Remove use of tensor library in B assembly
Martin Almquist <malmquist@stanford.edu>
parents: 963
diff changeset
376
1118
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
377 % Direction 1
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
378 for k = 1:m(1)
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
379 c = sparse(m(1),1);
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
380 c(k) = 1;
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
381 [~, B_1D] = ops{1}.D2(c);
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
382 for l = 1:m(2)
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
383 p = ind(:,l);
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
384 B{1}{(k-1)*m(2) + l}(p, p) = B_1D;
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
385 end
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
386 end
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
387
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
388 % Direction 2
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
389 for k = 1:m(2)
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
390 c = sparse(m(2),1);
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
391 c(k) = 1;
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
392 [~, B_1D] = ops{2}.D2(c);
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
393 for l = 1:m(1)
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
394 p = ind(l,:);
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
395 B{2}{(l-1)*m(2) + k}(p, p) = B_1D;
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
396 end
861
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
397 end
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
398 end
1118
07d0caf915e4 Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
Martin Almquist <malmquist@stanford.edu>
parents: 1060
diff changeset
399 obj.B = B;
965
db3411264b96 Remove use of tensor library in B assembly
Martin Almquist <malmquist@stanford.edu>
parents: 963
diff changeset
400
db3411264b96 Remove use of tensor library in B assembly
Martin Almquist <malmquist@stanford.edu>
parents: 963
diff changeset
401
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
402 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
403
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
404
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
405 % 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
406 % 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
407 % 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
408 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
1056
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
409 % on the first component. Can also be e.g.
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
410 % {'normal', 'd'} or {'tangential', 't'} for conditions on
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
411 % tangential/normal component.
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
412 % 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
413 % neighbour_scheme is an instance of Scheme that should be interfaced to.
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
414 % 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
415 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
416 default_arg('tuning', 1.2);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
417
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
418 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
419 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
420 type = bc{2};
1056
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
421 if ischar(comp)
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
422 comp = obj.getComponent(comp, boundary);
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
423 end
729
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
424
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
425 e = obj.getBoundaryOperatorForScalarField('e', boundary);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
426 tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
427 T = obj.getBoundaryTractionOperator(boundary);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
428 alpha = obj.getBoundaryOperatorForScalarField('alpha', boundary);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
429 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
430
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
431 E = obj.E;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
432 Hi = obj.Hi;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
433 LAMBDA = obj.LAMBDA;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
434 MU = obj.MU;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
435 RHOi = obj.RHOi;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
436
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
437 dim = obj.dim;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
438 m_tot = obj.grid.N();
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
439
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
440 % Preallocate
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
441 [~, col] = size(tau);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
442 closure = sparse(dim*m_tot, dim*m_tot);
1058
84933722ec0e Remove superfluous method get_boundary_number in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1057
diff changeset
443 penalty = sparse(dim*m_tot, col);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
444
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
445 k = comp;
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
446 switch type
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
447
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
448 % 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
449 case {'D','d','dirichlet','Dirichlet'}
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
450
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
451 % 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
452 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
453 C = transpose(T{k,i});
974
1c334842bf23 Extract tuning from alpha.
Martin Almquist <malmquist@stanford.edu>
parents: 973
diff changeset
454 A = -tuning*e*transpose(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
455 B = A + e*C;
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
456 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
457 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
458 end
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
459
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
460 % 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
461 case {'F','f','Free','free','traction','Traction','t','T'}
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
462 closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau';
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
463 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
464
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
465 % 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
466 otherwise
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
467 error('No such boundary condition: type = %s',type);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
468 end
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 % type Struct that specifies the interface coupling.
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
472 % Fields:
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
473 % -- tuning: penalty strength, defaults to 1.2
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
474 % -- interpolation: type of interpolation, default 'none'
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
475 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
476
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
477 defaultType.tuning = 1.2;
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
478 defaultType.interpolation = 'none';
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
479 default_struct('type', defaultType);
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
480
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
481 switch type.interpolation
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
482 case {'none', ''}
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
483 [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
484 case {'op','OP'}
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
485 [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
486 otherwise
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
487 error('Unknown type of interpolation: %s ', type.interpolation);
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
488 end
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
489 end
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
490
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
491 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
963
c75ddd568fcc Turn alpha into a boundary operator. Add properties H_w etc for getBoundaryQuadrature to work.
Martin Almquist <malmquist@stanford.edu>
parents: 961
diff changeset
492 tuning = type.tuning;
c75ddd568fcc Turn alpha into a boundary operator. Add properties H_w etc for getBoundaryQuadrature to work.
Martin Almquist <malmquist@stanford.edu>
parents: 961
diff changeset
493
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
494 % u denotes the solution in the own domain
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
495 % v denotes the solution in the neighbour domain
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
496 % Operators without subscripts are from the own domain.
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
497
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
498 % Get boundary operators
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
499 e = obj.getBoundaryOperator('e', boundary);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
500 tau = obj.getBoundaryOperator('tau', boundary);
1057
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
501
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
502 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
503 tau_v = neighbour_scheme.getBoundaryOperator('tau', neighbour_boundary);
1057
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
504
972
104f0af001e0 Clean up Elastic2dVariable.interfaceStandard
Martin Almquist <malmquist@stanford.edu>
parents: 970
diff changeset
505 H_gamma = obj.getBoundaryQuadrature(boundary);
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
506
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
507 % Operators and quantities that correspond to the own domain only
972
104f0af001e0 Clean up Elastic2dVariable.interfaceStandard
Martin Almquist <malmquist@stanford.edu>
parents: 970
diff changeset
508 Hi = obj.Hi_kron;
104f0af001e0 Clean up Elastic2dVariable.interfaceStandard
Martin Almquist <malmquist@stanford.edu>
parents: 970
diff changeset
509 RHOi = obj.RHOi_kron;
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
510
972
104f0af001e0 Clean up Elastic2dVariable.interfaceStandard
Martin Almquist <malmquist@stanford.edu>
parents: 970
diff changeset
511 % Penalty strength operators
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
512 alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha', boundary);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
513 alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha', neighbour_boundary);
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
514
972
104f0af001e0 Clean up Elastic2dVariable.interfaceStandard
Martin Almquist <malmquist@stanford.edu>
parents: 970
diff changeset
515 closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e');
104f0af001e0 Clean up Elastic2dVariable.interfaceStandard
Martin Almquist <malmquist@stanford.edu>
parents: 970
diff changeset
516 penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v');
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
517
972
104f0af001e0 Clean up Elastic2dVariable.interfaceStandard
Martin Almquist <malmquist@stanford.edu>
parents: 970
diff changeset
518 closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau';
104f0af001e0 Clean up Elastic2dVariable.interfaceStandard
Martin Almquist <malmquist@stanford.edu>
parents: 970
diff changeset
519 penalty = penalty - 1/2*RHOi*Hi*e*H_gamma*tau_v';
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
520
972
104f0af001e0 Clean up Elastic2dVariable.interfaceStandard
Martin Almquist <malmquist@stanford.edu>
parents: 970
diff changeset
521 closure = closure + 1/2*RHOi*Hi*tau*H_gamma*e';
104f0af001e0 Clean up Elastic2dVariable.interfaceStandard
Martin Almquist <malmquist@stanford.edu>
parents: 970
diff changeset
522 penalty = penalty - 1/2*RHOi*Hi*tau*H_gamma*e_v';
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
523
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
524 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
525
960
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
526 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
527 error('Non-conforming interfaces not implemented yet.');
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
528 end
ac566f3dc9b3 Add type to Elastic2dVariable.interface
Martin Almquist <malmquist@stanford.edu>
parents: 959
diff changeset
529
1056
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
530 % Returns the component number that is the tangential/normal component
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
531 % at the specified boundary
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
532 function comp = getComponent(obj, comp_str, boundary)
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
533 assertIsMember(comp_str, {'normal', 'tangential'});
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
534 assertIsMember(boundary, {'w', 'e', 's', 'n'});
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
535
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
536 switch boundary
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
537 case {'w', 'e'}
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
538 switch comp_str
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
539 case 'normal'
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
540 comp = 1;
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
541 case 'tangential'
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
542 comp = 2;
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
543 end
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
544 case {'s', 'n'}
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
545 switch comp_str
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
546 case 'normal'
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
547 comp = 2;
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
548 case 'tangential'
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
549 comp = 1;
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
550 end
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
551 end
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
552 end
b4fa176b4287 Elastic2dVariable: boundary condition type can now be a string, normal or tangential, as well as component number.
Martin Almquist <malmquist@stanford.edu>
parents: 974
diff changeset
553
726
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
554 % Returns the boundary operator op for the boundary specified by the string boundary.
1057
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
555 % op -- string
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
556 function o = getBoundaryOperator(obj, op, boundary)
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
557 assertIsMember(boundary, {'w', 'e', 's', 'n'})
1060
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
558 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'alpha', 'alpha1', 'alpha2'})
726
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
559
1057
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
560 switch op
970
23d9ca6755be Add getBoundaryQuadrature in Elastic2dVariable. Rename get_boundary_operator -> getBoundaryOperator. Add operators in getBoundaryOperator, with full size so that they work with multiblock.DiffOp.getBoundaryOperator.
Martin Almquist <malmquist@stanford.edu>
parents: 968
diff changeset
561
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
562 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
563 o = obj.([op, '_', boundary]);
861
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
564
1060
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
565 % Yields vector-valued penalty strength given displacement BC on all components
1057
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
566 case 'alpha'
1060
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
567 e = obj.getBoundaryOperator('e', boundary);
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
568 e_scalar = obj.getBoundaryOperatorForScalarField('e', boundary);
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
569 alpha_scalar = obj.getBoundaryOperatorForScalarField('alpha', boundary);
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
570 E = obj.E;
1060
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
571 [m, n] = size(alpha_scalar{1,1});
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
572 alpha = sparse(m*obj.dim, n*obj.dim);
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
573 for i = 1:obj.dim
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
574 for l = 1:obj.dim
1060
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
575 alpha = alpha + (e'*E{i}*e_scalar*alpha_scalar{i,l}'*E{l}')';
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
576 end
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
577 end
1060
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
578 o = alpha;
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
579
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
580 % Yields penalty strength for component 1 given displacement BC on all components
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
581 case 'alpha1'
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
582 alpha = obj.getBoundaryOperator('alpha', boundary);
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
583 e = obj.getBoundaryOperator('e', boundary);
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
584 e1 = obj.getBoundaryOperator('e1', boundary);
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
585
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
586 alpha1 = (e1'*e*alpha')';
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
587 o = alpha1;
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
588
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
589 % Yields penalty strength for component 2 given displacement BC on all components
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
590 case 'alpha2'
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
591 alpha = obj.getBoundaryOperator('alpha', boundary);
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
592 e = obj.getBoundaryOperator('e', boundary);
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
593 e2 = obj.getBoundaryOperator('e2', boundary);
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
594
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
595 alpha2 = (e2'*e*alpha')';
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
596 o = alpha2;
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
597 end
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
598
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
599 end
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
600
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
601 % Returns the boundary operator op for the boundary specified by the string boundary.
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
602 % op -- string
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
603 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
604 assertIsMember(boundary, {'w', 'e', 's', 'n'})
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
605 assertIsMember(op, {'e', 'alpha'})
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
606
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
607 switch op
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
608
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
609 case 'e'
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
610 o = obj.(['e_scalar', '_', boundary]);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
611
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
612 case 'alpha'
1060
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
613
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
614 % alpha{i,j} is the penalty strength on component i due to
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
615 % displacement BC for component j.
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
616 e = obj.getBoundaryOperatorForScalarField('e', boundary);
1057
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
617
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
618 LAMBDA = obj.LAMBDA;
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
619 MU = obj.MU;
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
620 dim = obj.dim;
861
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
621
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
622 switch boundary
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
623 case {'w', 'e'}
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
624 k = 1;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
625 case {'s', 'n'}
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
626 k = 2;
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
627 end
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
628
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
629 theta_R = obj.theta_R{k};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
630 theta_H = obj.theta_H{k};
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
631 theta_M = obj.theta_M{k};
861
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
Martin Almquist <malmquist@stanford.edu>
parents: 855
diff changeset
632
1057
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
633 a_lambda = dim/theta_H + 1/theta_R;
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
634 a_mu_i = 2/theta_M;
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
635 a_mu_ij = 2/theta_H + 1/theta_R;
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
636
1057
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
637 d = @kroneckerDelta; % Kronecker delta
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
638 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
639
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
640 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
641 + d(i,j)* a_mu_i*MU ...
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
642 + db(i,j)*a_mu_ij*MU;
1060
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
643
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
644 alpha = cell(obj.dim, obj.dim);
1057
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
645 for i = 1:obj.dim
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
646 for j = 1:obj.dim
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
647 alpha{i,j} = d(i,j)*alpha_func(i,k)*e;
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
648 end
1057
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
649 end
ff274c7404cc In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
Martin Almquist <malmquist@stanford.edu>
parents: 1056
diff changeset
650 o = alpha;
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
651 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
652
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
653 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
654
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
655 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
1060
e40899094f20 Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
Martin Almquist <malmquist@stanford.edu>
parents: 1059
diff changeset
656 % Formula: tau_i = T_ij u_j
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
657 % op -- string
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
658 function T = getBoundaryTractionOperator(obj, boundary)
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
659 assertIsMember(boundary, {'w', 'e', 's', 'n'})
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
660
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
661 T = obj.(['T', '_', boundary]);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
662 end
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
663
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
664 % Returns square boundary quadrature matrix, of dimension
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
665 % corresponding to the number of boundary unknowns
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
666 %
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
667 % boundary -- string
970
23d9ca6755be Add getBoundaryQuadrature in Elastic2dVariable. Rename get_boundary_operator -> getBoundaryOperator. Add operators in getBoundaryOperator, with full size so that they work with multiblock.DiffOp.getBoundaryOperator.
Martin Almquist <malmquist@stanford.edu>
parents: 968
diff changeset
668 function H = getBoundaryQuadrature(obj, boundary)
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
669 assertIsMember(boundary, {'w', 'e', 's', 'n'})
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
670
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
671 H = obj.getBoundaryQuadratureForScalarField(boundary);
970
23d9ca6755be Add getBoundaryQuadrature in Elastic2dVariable. Rename get_boundary_operator -> getBoundaryOperator. Add operators in getBoundaryOperator, with full size so that they work with multiblock.DiffOp.getBoundaryOperator.
Martin Almquist <malmquist@stanford.edu>
parents: 968
diff changeset
672 I_dim = speye(obj.dim, obj.dim);
23d9ca6755be Add getBoundaryQuadrature in Elastic2dVariable. Rename get_boundary_operator -> getBoundaryOperator. Add operators in getBoundaryOperator, with full size so that they work with multiblock.DiffOp.getBoundaryOperator.
Martin Almquist <malmquist@stanford.edu>
parents: 968
diff changeset
673 H = kron(H, I_dim);
23d9ca6755be Add getBoundaryQuadrature in Elastic2dVariable. Rename get_boundary_operator -> getBoundaryOperator. Add operators in getBoundaryOperator, with full size so that they work with multiblock.DiffOp.getBoundaryOperator.
Martin Almquist <malmquist@stanford.edu>
parents: 968
diff changeset
674 end
23d9ca6755be Add getBoundaryQuadrature in Elastic2dVariable. Rename get_boundary_operator -> getBoundaryOperator. Add operators in getBoundaryOperator, with full size so that they work with multiblock.DiffOp.getBoundaryOperator.
Martin Almquist <malmquist@stanford.edu>
parents: 968
diff changeset
675
1059
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
676 % Returns square boundary quadrature matrix, of dimension
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
677 % corresponding to the number of boundary grid points
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
678 %
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
679 % boundary -- string
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
680 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
681 assertIsMember(boundary, {'w', 'e', 's', 'n'})
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
682
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
683 H_b = obj.(['H_', boundary]);
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
684 end
d8ab528f10f6 Massive boundary operator cleanup in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 1058
diff changeset
685
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
686 function N = size(obj)
734
eebe24a636c7 Make Elastic2dVariable.size account for components.
Martin Almquist <malmquist@stanford.edu>
parents: 730
diff changeset
687 N = obj.dim*prod(obj.m);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
688 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
689 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
690 end