Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dVariable.m @ 1031:2ef20d00b386 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
For easier comparison, return both the first order and residual viscosity when evaluating the residual. Add the first order and residual viscosity to the state of the RungekuttaRV time steppers
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 17 Jan 2019 10:25:06 +0100
parents 21394c78c72e
children a4ad90b37998 a2fcc4cf2298
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
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 % Diagonal matrices for varible coefficients
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18 LAMBDA % Variable coefficient, related to dilation
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 MU % Shear modulus, variable coefficient
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 RHO, RHOi % Density, variable
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 D % Total operator
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 D1 % First derivatives
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25 % Second derivatives
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 D2_lambda
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 D2_mu
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29 % Traction operators used for BC
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 T_l, T_r
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31 tau_l, tau_r
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 H, Hi % Inner products
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
34
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35 phi % Borrowing constant for (d1 - e^T*D1) from R
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36 gamma % Borrowing constant for d1 from M
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 H11 % First element of H
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
38
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
39 % Borrowing from H, M, and R
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
40 thH
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
41 thM
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
42 thR
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
43
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 e_l, e_r
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45 d1_l, d1_r % Normal derivatives at the boundary
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 E % E{i}^T picks out component i
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
47
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 H_boundary % Boundary inner products
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 % Kroneckered norms and coefficients
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 RHOi_kron
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 Hi_kron
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55 methods
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet)
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 default_arg('lambda_fun', @(x,y) 0*x+1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 default_arg('mu_fun', @(x,y) 0*x+1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 default_arg('rho_fun', @(x,y) 0*x+1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 dim = 2;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 assert(isa(g, 'grid.Cartesian'))
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 lambda = grid.evalOn(g, lambda_fun);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 mu = grid.evalOn(g, mu_fun);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 rho = grid.evalOn(g, rho_fun);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 m = g.size();
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70 m_tot = g.N();
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 h = g.scaling();
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 lim = g.lim;
729
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
74 if isempty(lim)
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
75 x = g.x;
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
76 lim = cell(length(x),1);
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
77 for i = 1:length(x)
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
78 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
79 end
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
80 end
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 % 1D operators
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 ops = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 ops{i} = opSet{i}(m(i), lim{i}, order);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 end
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 % Borrowing constants
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 beta = ops{i}.borrowing.R.delta_D;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 obj.H11{i} = ops{i}.borrowing.H11;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 obj.phi{i} = beta/obj.H11{i};
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 obj.gamma{i} = ops{i}.borrowing.M.d1;
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
94
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
95 % Better names
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
96 obj.thR{i} = ops{i}.borrowing.R.delta_D;
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
97 obj.thM{i} = ops{i}.borrowing.M.d1;
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
98 obj.thH{i} = ops{i}.borrowing.H11;
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101 I = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 D1 = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 D2 = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 H = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 Hi = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 e_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 e_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 d1_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 d1_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 I{i} = speye(m(i));
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113 D1{i} = ops{i}.D1;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 D2{i} = ops{i}.D2;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 H{i} = ops{i}.H;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 Hi{i} = ops{i}.HI;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 e_l{i} = ops{i}.e_l;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118 e_r{i} = ops{i}.e_r;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119 d1_l{i} = ops{i}.d1_l;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 d1_r{i} = ops{i}.d1_r;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 %====== Assemble full operators ========
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 LAMBDA = spdiag(lambda);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 obj.LAMBDA = LAMBDA;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 MU = spdiag(mu);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 obj.MU = MU;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128 RHO = spdiag(rho);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129 obj.RHO = RHO;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 obj.RHOi = inv(RHO);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132 obj.D1 = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 obj.D2_lambda = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134 obj.D2_mu = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135 obj.e_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 obj.e_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 obj.d1_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 obj.d1_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140 % D1
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141 obj.D1{1} = kron(D1{1},I{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142 obj.D1{2} = kron(I{1},D1{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144 % Boundary operators
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145 obj.e_l{1} = kron(e_l{1},I{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 obj.e_l{2} = kron(I{1},e_l{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 obj.e_r{1} = kron(e_r{1},I{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148 obj.e_r{2} = kron(I{1},e_r{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
150 obj.d1_l{1} = kron(d1_l{1},I{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
151 obj.d1_l{2} = kron(I{1},d1_l{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
152 obj.d1_r{1} = kron(d1_r{1},I{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 obj.d1_r{2} = kron(I{1},d1_r{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155 % D2
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
156 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
157 obj.D2_lambda{i} = sparse(m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
158 obj.D2_mu{i} = sparse(m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
160 ind = grid.funcToMatrix(g, 1:m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
161
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
162 for i = 1:m(2)
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
163 D_lambda = D2{1}(lambda(ind(:,i)));
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
164 D_mu = D2{1}(mu(ind(:,i)));
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
165
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
166 p = ind(:,i);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
167 obj.D2_lambda{1}(p,p) = D_lambda;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168 obj.D2_mu{1}(p,p) = D_mu;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
169 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
170
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
171 for i = 1:m(1)
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
172 D_lambda = D2{2}(lambda(ind(i,:)));
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173 D_mu = D2{2}(mu(ind(i,:)));
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 p = ind(i,:);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
176 obj.D2_lambda{2}(p,p) = D_lambda;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
177 obj.D2_mu{2}(p,p) = D_mu;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 % Quadratures
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181 obj.H = kron(H{1},H{2});
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
182 obj.Hi = inv(obj.H);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
183 obj.H_boundary = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184 obj.H_boundary{1} = H{2};
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
185 obj.H_boundary{2} = H{1};
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187 % E{i}^T picks out component i.
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
188 E = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
189 I = speye(m_tot,m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191 e = sparse(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192 e(i) = 1;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
193 E{i} = kron(I,e);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
194 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195 obj.E = E;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
197 % Differentiation matrix D (without SAT)
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198 D2_lambda = obj.D2_lambda;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 D2_mu = obj.D2_mu;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200 D1 = obj.D1;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201 D = sparse(dim*m_tot,dim*m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 d = @kroneckerDelta; % Kronecker delta
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 for j = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206 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
207 db(i,j)*D1{i}*LAMBDA*D1{j}*E{j}' ...
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 );
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 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
210 db(i,j)*D1{j}*MU*D1{i}*E{j}' + ...
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 D2_mu{j}*E{i}' ...
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212 );
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 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 obj.D = D;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 %=========================================%
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
218 % Numerical traction operators for BC.
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219 % Because d1 =/= e0^T*D1, the numerical tractions are different
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 % at every boundary.
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221 T_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 T_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 tau_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 tau_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 % tau^{j}_i = sum_k T^{j}_{ik} u_k
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
227 d1_l = obj.d1_l;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228 d1_r = obj.d1_r;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 e_l = obj.e_l;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 e_r = obj.e_r;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231 D1 = obj.D1;
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 % Loop over boundaries
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
234 for j = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
235 T_l{j} = cell(dim,dim);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
236 T_r{j} = cell(dim,dim);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
237 tau_l{j} = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
238 tau_r{j} = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
239
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 % Loop over components
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
242 tau_l{j}{i} = sparse(m_tot,dim*m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
243 tau_r{j}{i} = sparse(m_tot,dim*m_tot);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
244 for k = 1:dim
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
245 T_l{j}{i,k} = ...
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
246 -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})...
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
247 -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})...
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248 -d(i,k)*MU*e_l{j}*d1_l{j}';
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
249
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
250 T_r{j}{i,k} = ...
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
251 d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})...
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
252 +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})...
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253 +d(i,k)*MU*e_r{j}*d1_r{j}';
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
254
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
257 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
258
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
259 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
260 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
261 obj.T_l = T_l;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
262 obj.T_r = T_r;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
263 obj.tau_l = tau_l;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
264 obj.tau_r = tau_r;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
265
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
266 % Kroneckered norms and coefficients
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
267 I_dim = speye(dim);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
268 obj.RHOi_kron = kron(obj.RHOi, I_dim);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
269 obj.Hi_kron = kron(obj.Hi, I_dim);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
271 % Misc.
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
272 obj.m = m;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
273 obj.h = h;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
274 obj.order = order;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
275 obj.grid = g;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
276 obj.dim = dim;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
277
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
278 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
279
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
280
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
281 % 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
282 % 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
283 % 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
284 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
285 % on the first component.
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
286 % 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
287 % neighbour_scheme is an instance of Scheme that should be interfaced to.
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
288 % 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
289 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
290 default_arg('tuning', 1.2);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
291
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
292 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
293 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
294 type = bc{2};
729
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
295
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
296 % j is the coordinate direction of the boundary
726
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
297 j = obj.get_boundary_number(boundary);
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
298 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
299
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
300 E = obj.E;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
301 Hi = obj.Hi;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
302 LAMBDA = obj.LAMBDA;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
303 MU = obj.MU;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
304 RHOi = obj.RHOi;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
305
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
306 dim = obj.dim;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
307 m_tot = obj.grid.N();
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
308
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
309 % Preallocate
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
310 closure = sparse(dim*m_tot, dim*m_tot);
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
311 penalty = sparse(dim*m_tot, m_tot/obj.m(j));
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
312
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
313 k = comp;
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
314 switch type
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
315
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
316 % 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
317 case {'D','d','dirichlet','Dirichlet'}
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
318
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
319 phi = obj.phi{j};
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
320 h = obj.h(j);
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
321 h11 = obj.H11{j}*h;
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
322 gamma = obj.gamma{j};
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
323
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
324 a_lambda = dim/h11 + 1/(h11*phi);
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
325 a_mu_i = 2/(gamma*h);
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
326 a_mu_ij = 2/h11 + 1/(h11*phi);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
327
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
328 d = @kroneckerDelta; % Kronecker delta
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
329 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
330 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
331 + d(i,j)* a_mu_i*MU ...
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
332 + db(i,j)*a_mu_ij*MU );
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
333
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
334 % 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
335 for i = 1:dim
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
336 C = T{k,i};
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
337 A = -d(i,k)*alpha(i,j);
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
338 B = A + C;
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
339 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
340 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
341 end
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
342
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
343 % 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
344 case {'F','f','Free','free','traction','Traction','t','T'}
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
345 closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} );
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
346 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
347
795
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
348 % 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
349 otherwise
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible.
Martin Almquist <malmquist@stanford.edu>
parents: 738
diff changeset
350 error('No such boundary condition: type = %s',type);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
353
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
354 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
355 % u denotes the solution in the own domain
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
356 % v denotes the solution in the neighbour domain
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
357 % Operators without subscripts are from the own domain.
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
358 tuning = 1.2;
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
359
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
360 % j is the coordinate direction of the boundary
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
361 j = obj.get_boundary_number(boundary);
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
362 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
363
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
364 % Get boundary operators
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
365 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
366 [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
367
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
368 % Operators and quantities that correspond to the own domain only
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
369 Hi = obj.Hi;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
370 RHOi = obj.RHOi;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
371 dim = obj.dim;
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
372
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
373 %--- Other operators ----
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
374 m_tot_u = obj.grid.N();
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
375 E = obj.E;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
376 LAMBDA_u = obj.LAMBDA;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
377 MU_u = obj.MU;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
378 lambda_u = e'*LAMBDA_u*e;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
379 mu_u = e'*MU_u*e;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
380
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
381 m_tot_v = neighbour_scheme.grid.N();
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
382 E_v = neighbour_scheme.E;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
383 LAMBDA_v = neighbour_scheme.LAMBDA;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
384 MU_v = neighbour_scheme.MU;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
385 lambda_v = e_v'*LAMBDA_v*e_v;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
386 mu_v = e_v'*MU_v*e_v;
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
387 %-------------------------
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
388
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
389 % Borrowing constants
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
390 h_u = obj.h(j);
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
391 thR_u = obj.thR{j}*h_u;
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
392 thM_u = obj.thM{j}*h_u;
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
393 thH_u = obj.thH{j}*h_u;
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
394
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
395 h_v = neighbour_scheme.h(j_v);
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
396 thR_v = neighbour_scheme.thR{j_v}*h_v;
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
397 thH_v = neighbour_scheme.thH{j_v}*h_v;
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
398 thM_v = neighbour_scheme.thM{j_v}*h_v;
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
399
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
400 % alpha = penalty strength for normal component, beta for tangential
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
401 alpha_u = dim*lambda_u/(4*thH_u) + lambda_u/(4*thR_u) + mu_u/(2*thM_u);
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
402 alpha_v = dim*lambda_v/(4*thH_v) + lambda_v/(4*thR_v) + mu_v/(2*thM_v);
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
403 beta_u = mu_u/(2*thH_u) + mu_u/(4*thR_u);
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
404 beta_v = mu_v/(2*thH_v) + mu_v/(4*thR_v);
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
405 alpha = alpha_u + alpha_v;
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
406 beta = beta_u + beta_v;
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
407
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
408 d = @kroneckerDelta; % Kronecker delta
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
409 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
410 strength = @(i,j) tuning*(d(i,j)*alpha + db(i,j)*beta);
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
411
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
412 % Preallocate
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
413 closure = sparse(dim*m_tot_u, dim*m_tot_u);
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
414 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
415
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
416 % Loop over components that penalties end up on
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
417 for i = 1:dim
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
418 closure = closure - E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e'*E{i}';
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
419 penalty = penalty + E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e_v'*E_v{i}';
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
420
729
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
421 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
aa8cf3851de8 Update multiblock.DiffOp to work for systems.
Martin Almquist <malmquist@stanford.edu>
parents: 727
diff changeset
422 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
423
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
424 % Loop over components that we have interface conditions on
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
425 for k = 1:dim
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
426 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}';
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
427 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}';
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
428 end
727
6d5953fc090e First implementation of elastic interface
Martin Almquist <malmquist@stanford.edu>
parents: 726
diff changeset
429 end
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
430 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
431
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
432 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
433 function [j, nj] = get_boundary_number(obj, boundary)
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
434
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
435 switch boundary
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
436 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
437 j = 1;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
438 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
439 j = 2;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
440 otherwise
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
441 error('No such boundary: boundary = %s',boundary);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
442 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
443
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
444 switch boundary
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
445 case {'w','W','west','West','s','S','south','South'}
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
446 nj = -1;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
447 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
448 nj = 1;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
449 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
450 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
451
726
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
452 % Returns the boundary operator op for the boundary specified by the string boundary.
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
453 % op: may be a cell array of strings
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
454 function [varargout] = get_boundary_operator(obj, op, boundary)
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
455
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
456 switch boundary
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
457 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
458 j = 1;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
459 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
460 j = 2;
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
461 otherwise
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
462 error('No such boundary: boundary = %s',boundary);
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
463 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
464
726
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
465 if ~iscell(op)
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
466 op = {op};
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
467 end
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
468
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
469 for i = 1:length(op)
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
470 switch op{i}
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
471 case 'e'
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
472 switch boundary
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
473 case {'w','W','west','West','s','S','south','South'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
474 varargout{i} = obj.e_l{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
475 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
476 varargout{i} = obj.e_r{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
477 end
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
478 case 'd'
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
479 switch boundary
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
480 case {'w','W','west','West','s','S','south','South'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
481 varargout{i} = obj.d1_l{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
482 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
483 varargout{i} = obj.d1_r{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
484 end
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
485 case 'H'
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
486 varargout{i} = obj.H_boundary{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
487 case 'T'
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
488 switch boundary
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
489 case {'w','W','west','West','s','S','south','South'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
490 varargout{i} = obj.T_l{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
491 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
492 varargout{i} = obj.T_r{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
493 end
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
494 case 'tau'
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
495 switch boundary
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
496 case {'w','W','west','West','s','S','south','South'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
497 varargout{i} = obj.tau_l{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
498 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
499 varargout{i} = obj.tau_r{j};
813
b374a8aa9246 Correct interface penalty strength in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 795
diff changeset
500 end
726
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
501 otherwise
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
502 error(['No such operator: operator = ' op{i}]);
37d5d69b1a0d Refactor BC code in Elastic2dVariable
Martin Almquist <malmquist@stanford.edu>
parents: 689
diff changeset
503 end
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
504 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
505
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
506 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
507
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
508 function N = size(obj)
734
eebe24a636c7 Make Elastic2dVariable.size account for components.
Martin Almquist <malmquist@stanford.edu>
parents: 730
diff changeset
509 N = obj.dim*prod(obj.m);
687
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
510 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
511 end
e8fc3aa1faf6 Rename elastic scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
512 end