Mercurial > repos > public > sbplib

annotate +scheme/ViscoElastic2d.m @ 1313:a4c8bf2371f3 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add viscoelastic fault interface
author Martin Almquist <malmquist@stanford.edu>
date Wed, 22 Jul 2020 15:35:25 -0700
parents 6d64b57caf46
children 58df4a35fe43
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 classdef ViscoElastic2d < scheme.Scheme
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 % Discretizes the visco-elastic wave equation in curvilinear coordinates.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 % Assumes fully compatible operators.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 properties
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 m % Number of points in each direction, possibly a vector
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 h % Grid spacing
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 grid
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13 order % Order of accuracy for the approximation
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 % Diagonal matrices for variable coefficients
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16 % J, Ji
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 RHO % Density
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18 C % Elastic stiffness tensor
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 ETA % Effective viscosity, used in strain rate eq
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21 D % Total operator
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 Delastic % Elastic operator (momentum balance)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 Dviscous % Acts on viscous strains in momentum balance
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24 DstrainRate % Acts on u and gamma, returns strain rate gamma_t
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 D1, D1Tilde % Physical derivatives
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 sigma % Cell matrix of physical stress operators
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29 % Inner products
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 H
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32 % Boundary inner products (for scalar field)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 H_w, H_e, H_s, H_n
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35 % Restriction operators
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36 Eu, Egamma % Pick out all components of u/gamma
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 eU, eGamma % Pick out one specific component
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39 % Bundary restriction ops
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 n_w, n_e, n_s, n_n % Physical normals
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
43 tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45 elasticObj
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 methods
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 % The coefficients can either be function handles or grid functions
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 function obj = ViscoElastic2d(g, order, rho, C, eta)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 default_arg('rho', @(x,y) 0*x+1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 default_arg('eta', @(x,y) 0*x+1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55 dim = 2;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 C_default = cell(dim,dim,dim,dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 C_default{i,j,k,l} = @(x,y) 0*x ;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 default_arg('C', C_default);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 assert(isa(g, 'grid.Curvilinear'));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 if isa(rho, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 rho = grid.evalOn(g, rho);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 if isa(eta, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 eta = grid.evalOn(g, eta);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 C_mat = cell(dim,dim,dim,dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 if isa(C{i,j,k,l}, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 obj.C = C_mat;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 elasticObj = scheme.Elastic2dCurvilinearAnisotropic(g, order, rho, C);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 % Construct a pair of first derivatives
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 K = elasticObj.K;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 K{i,j} = spdiag(K{i,j});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 J = elasticObj.J;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 Ji = elasticObj.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 D_ref = elasticObj.refObj.D1;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 D1 = cell(dim, 1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 D1Tilde = cell(dim, 1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 D1{i} = 0*D_ref{i};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 D1Tilde{i} = 0*D_ref{i};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 D1{i} = D1{i} + K{i,j}*D_ref{j};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113 D1Tilde{i} = D1Tilde{i} + Ji*D_ref{j}*J*K{i,j};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 obj.D1 = D1;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 obj.D1Tilde = D1Tilde;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119 eU = elasticObj.E;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 % Storage order for gamma: 11-12-21-22.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 I = speye(g.N(), g.N());
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 eGamma = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 e = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 e{1,1} = [1;0;0;0];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 e{1,2} = [0;1;0;0];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 e{2,1} = [0;0;1;0];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128 e{2,2} = [0;0;0;1];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 eGamma{i,j} = kron(I, e{i,j});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135 % Store u first, then gamma
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 mU = dim*g.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 mGamma = dim^2*g.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 Iu = speye(mU, mU);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139 Igamma = speye(mGamma, mGamma);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141 Eu = cell2mat({Iu, sparse(mU, mGamma)})';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142 Egamma = cell2mat({sparse(mGamma, mU), Igamma})';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145 eU{i} = Eu*eU{i};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149 eGamma{i,j} = Egamma*eGamma{i,j};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
150 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
151 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
152
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 obj.eGamma = eGamma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154 obj.eU = eU;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155 obj.Egamma = Egamma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
156 obj.Eu = Eu;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
157
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
158 % Build stress operator
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 sigma = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
160 C = obj.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
161 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
162 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
163 sigma{i,j} = spalloc(g.N(), (dim^2 + dim)*g.N(), order^2*g.N());
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
164 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
165 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
166 sigma{i,j} = sigma{i,j} + C{i,j,k,l}*(D1{k}*eU{l}' - eGamma{k,l}');
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
167 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
169 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
170 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
171
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
172 % Elastic operator
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173 Delastic = Eu*elasticObj.D*Eu';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 % Add viscous strains to momentum balance
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
176 RHOi = spdiag(1./rho);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
177 Dviscous = spalloc((dim^2 + dim)*g.N(), (dim^2 + dim)*g.N(), order^2*(dim^2 + dim)*g.N());
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
182 Dviscous = Dviscous - eU{j}*RHOi*D1Tilde{i}*C{i,j,k,l}*eGamma{k,l}';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
183 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
185 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
188 ETA = spdiag(eta);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
189 DstrainRate = 0*Delastic;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192 DstrainRate = DstrainRate + eGamma{i,j}*(ETA\sigma{i,j});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
193 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
194 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196 obj.D = Delastic + Dviscous + DstrainRate;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
197 obj.Delastic = Delastic;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198 obj.Dviscous = Dviscous;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 obj.DstrainRate = DstrainRate;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200 obj.sigma = sigma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 %---- Set remaining object properties ------
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 obj.RHO = elasticObj.RHO;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 obj.ETA = ETA;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 obj.H = elasticObj.H;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207 obj.n_w = elasticObj.n_w;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 obj.n_e = elasticObj.n_e;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 obj.n_s = elasticObj.n_s;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210 obj.n_n = elasticObj.n_n;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
212 obj.tangent_w = elasticObj.tangent_w;
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
213 obj.tangent_e = elasticObj.tangent_e;
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
214 obj.tangent_s = elasticObj.tangent_s;
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
215 obj.tangent_n = elasticObj.tangent_n;
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
216
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 obj.H_w = elasticObj.H_w;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
218 obj.H_e = elasticObj.H_e;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219 obj.H_s = elasticObj.H_s;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 obj.H_n = elasticObj.H_n;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 obj.e_scalar_w = elasticObj.e_scalar_w;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 obj.e_scalar_e = elasticObj.e_scalar_e;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 obj.e_scalar_s = elasticObj.e_scalar_s;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 obj.e_scalar_n = elasticObj.e_scalar_n;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
227 % Misc.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228 obj.elasticObj = elasticObj;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 obj.m = elasticObj.m;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 obj.h = elasticObj.h;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
232 obj.order = order;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
233 obj.grid = g;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
234 obj.dim = dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
235
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
236 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
237
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
238
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
239 % Closure functions return the operators applied to the own domain to close the boundary
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
242 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
243 % on the first component. Can also be e.g.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
244 % {'normal', 'd'} or {'tangential', 't'} for conditions on
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
245 % tangential/normal component.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
246 % data is a function returning the data that should be applied at the boundary.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
247 % neighbour_scheme is an instance of Scheme that should be interfaced to.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248 % neighbour_boundary is a string specifying which boundary to interface to.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
249
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
250 % For displacement bc:
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
251 % bc = {comp, 'd', dComps},
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
252 % where
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253 % dComps = vector of components with displacement BC. Default: 1:dim.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
254 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256 default_arg('tuning', 1.0);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
257 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
258
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
259 component = bc{1};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
260 type = bc{2};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
261 dim = obj.dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
262
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
263 n = obj.getNormal(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
264 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
265 e = obj.getBoundaryOperatorForScalarField('e', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
266
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
267 H = obj.H;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
268 RHO = obj.RHO;
1309
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
269 ETA = obj.ETA;
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270 C = obj.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
271 Eu = obj.Eu;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
272 eU = obj.eU;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
273 eGamma = obj.eGamma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
274
1311
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
275 % Get elastic closure and penalty
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
276 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning);
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
277 closure = Eu*closure*Eu';
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
278 penalty = Eu*penalty;
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
279
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
280 switch component
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
281 case 't'
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
282 dir = obj.getTangent(boundary);
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
283 case 'n'
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
284 dir = obj.getNormal(boundary);
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
285 case 1
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
286 dir = {1, 0};
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
287 case 2
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
288 dir = {0, 1};
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
289 end
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
290
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
291 switch type
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
292 case {'F','f','Free','free','traction','Traction','t','T'}
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
293
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
294 % Add viscous part of closure
1310
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
295 for i = 1:dim
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
296 for j = 1:dim
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
297 for k = 1:dim
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
298 for l = 1:dim
1310
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
299 for m = 1:dim
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
300 closure = closure + eU{m}*( (RHO*H)\(C{i,j,k,l}*e*dir{m}*dir{j}*H_gamma*n{i}*e'*eGamma{k,l}') );
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
301 end
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
302 end
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
303 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
304 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
305 end
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
306
1309
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
307 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
308
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
309 % Add penalty to strain rate eq
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
310 for i = 1:dim
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
311 for j = 1:dim
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
312 for k = 1:dim
1310
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
313 for l = 1:dim
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
314 for m = 1:dim
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
315 closure = closure - eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*dir{l}*dir{m}*n{k}*e'*eU{m}') );
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
316 end
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
317 penalty = penalty + eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*dir{l}*n{k}) );
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
318 end
1309
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
319 end
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
320 end
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
321 end
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
322
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
323 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
324
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
325 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
326
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
327 % type Struct that specifies the interface coupling.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
328 % Fields:
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
329 % -- tuning: penalty strength, defaults to 1.0
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
330 % -- interpolation: type of interpolation, default 'none'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
331 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
332
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
333 defaultType.tuning = 1.0;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
334 defaultType.interpolation = 'none';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
335 defaultType.type = 'standard';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
336 default_struct('type', defaultType);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
337
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
338 switch type.type
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
339 case 'standard'
1312
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
340 [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
341 case 'frictionalFault'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
342 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
343 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
344
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
345 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
346
1312
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
347 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
348
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
349 % u denotes the solution in the own domain
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
350 % v denotes the solution in the neighbour domain
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
351 u = obj;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
352 v = neighbour_scheme;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
353
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
354 dim = obj.dim;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
355
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
356 n = u.getNormal(boundary);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
357 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
358 e = u.getBoundaryOperatorForScalarField('e', boundary);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
359
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
360 ev = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
361
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
362 H = u.H;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
363 RHO = u.RHO;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
364 ETA = u.ETA;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
365 C = u.C;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
366 Eu = u.Eu;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
367 eU = u.eU;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
368 eGamma = u.eGamma;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
369
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
370 CV = v.C;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
371 Ev = v.Eu;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
372 eV = v.eU;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
373 eGammaV = v.eGamma;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
374 nV = v.getNormal(neighbour_boundary);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
375
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
376
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
377 % Get elastic closure and penalty
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
378 [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
379 closure = Eu*closure*Eu';
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
380 penalty = Eu*penalty*Ev';
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
381
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
382 % Add viscous part of traction coupling
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
383 for i = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
384 for j = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
385 for k = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
386 for l = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
387 closure = closure + 1/2*eU{j}*( (RHO*H)\(C{i,j,k,l}*e*H_gamma*n{i}*e'*eGamma{k,l}') );
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
388 penalty = penalty + 1/2*eU{j}*( (RHO*H)\(e*H_gamma*nV{i}*(ev'*CV{i,j,k,l}*ev)*ev'*eGammaV{k,l}') );
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
389 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
390 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
391 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
392 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
393
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
394 % Add penalty to strain rate eq
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
395 for i = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
396 for j = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
397 for k = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
398 for l = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
399 closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}*e'*eU{l}') );
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
400 penalty = penalty + 1/2*eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}*ev'*eV{l}') );
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
401 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
402 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
403 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
404 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
405
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
406
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
407 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
408
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
409 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
410 tuning = type.tuning;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
411
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
412 % u denotes the solution in the own domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
413 % v denotes the solution in the neighbour domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
414 u = obj;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
415 v = neighbour_scheme;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
416
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
417 dim = obj.dim;
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
418
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
419 n = u.getNormal(boundary);
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
420 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
421 e = u.getBoundaryOperatorForScalarField('e', boundary);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
422
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
423 ev = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
424
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
425 H = u.H;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
426 RHO = u.RHO;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
427 ETA = u.ETA;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
428 C = u.C;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
429 Eu = u.Eu;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
430 eU = u.eU;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
431 eGamma = u.eGamma;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
432 Egamma = u.Egamma;
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
433
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
434 CV = v.C;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
435 Ev = v.Eu;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
436 eV = v.eU;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
437 eGammaV = v.eGamma;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
438 nV = v.getNormal(neighbour_boundary);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
439
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
440 % Reduce stiffness tensors to boundary size
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
441 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
442 for j = 1:dim
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
443 for k = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
444 for l = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
445 C{i,j,k,l} = e'*C{i,j,k,l}*e;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
446 CV{i,j,k,l} = ev'*CV{i,j,k,l}*ev;
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
447 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
448 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
449 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
450 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
451
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
452 % Get elastic closure and penalty
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
453 [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
454 closure = Eu*closure*Eu';
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
455 penalty = Eu*penalty*Ev';
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
456
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
457 % ---- Tangential tractions are imposed just like traction BC ------
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
458 % We only need the viscous part
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
459 closure_tangential = obj.boundary_condition(boundary, {'t', 't'});
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
460 closure = closure + closure_tangential*Egamma*Egamma';
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
461
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
462
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
463 % ------ Coupling of normal component -----------
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
464 % Add viscous part of traction coupling
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
465 for i = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
466 for j = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
467 for k = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
468 for l = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
469 for m = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
470 closure = closure + 1/2*eU{m}*( (RHO*H)\(e*n{m}*H_gamma*n{j}*n{i}*C{i,j,k,l}*e'*eGamma{k,l}') );
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
471 penalty = penalty - 1/2*eU{m}*( (RHO*H)\(e*n{m}*H_gamma*nV{j}*nV{i}*CV{i,j,k,l}*ev'*eGammaV{k,l}') );
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
472 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
473 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
474 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
475 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
476 end
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
477
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
478 % Add penalty to strain rate eq
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
479 for i = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
480 for j = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
481 for k = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
482 for l = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
483 for m = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
484 closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma*n{m}*e'*eU{m}') );
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
485 penalty = penalty - 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma*nV{m}*ev'*eV{m}') );
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
486 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
487 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
488 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
489 end
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
490 end
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
491 %-------------------------------------------------
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
492
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
493 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
494
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
495 % Returns the outward unit normal vector for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
496 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
497 function n = getNormal(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
498 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
499
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
500 n = obj.(['n_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
501 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
502
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
503 % Returns the unit tangent vector for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
504 % t is a cell of diagonal matrices for each normal component, t{1} = t_1, t{2} = t_2.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
505 function t = getTangent(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
506 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
507
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
508 t = obj.(['tangent_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
509 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
510
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
511 % Returns the boundary operator op for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
512 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
513 function o = getBoundaryOperator(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
514 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
515 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
516
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
517 o = obj.elasticObj.([op, '_', boundary]);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
518
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
519 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
520
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
521 % Returns the boundary operator op for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
522 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
523 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
524 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
525 assertIsMember(op, {'e'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
526
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
527 switch op
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
528
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
529 case 'e'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
530 o = obj.(['e_scalar', '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
531 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
532
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
533 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
534
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
535 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
536 % Formula: tau_i = T_ij u_j
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
537 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
538 function T = getBoundaryTractionOperator(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
539 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
540
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
541 T = obj.(['T', '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
542 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
543
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
544 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
545 % corresponding to the number of boundary unknowns
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
546 %
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
547 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
548 function H = getBoundaryQuadrature(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
549 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
550
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
551 H = obj.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
552 I_dim = speye(obj.dim, obj.dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
553 H = kron(H, I_dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
554 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
555
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
556 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
557 % corresponding to the number of boundary grid points
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
558 %
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
559 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
560 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
561 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
562
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
563 H_b = obj.(['H_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
564 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
565
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
566 function N = size(obj)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
567 N = (obj.dim + obj.dim^2)*prod(obj.m);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
568 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
569 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
570 end