Mercurial > repos > public > sbplib

annotate +scheme/ViscoElastic2d.m @ 1307:fcca6ad8b102 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add diffOp for viscoElastic
author Martin Almquist <malmquist@stanford.edu>
date Sun, 19 Jul 2020 20:30:16 -0700
parents
children 5016f3f3badb
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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 elasticObj
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 methods
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 % The coefficients can either be function handles or grid functions
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 function obj = ViscoElastic2d(g, order, rho, C, eta)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 default_arg('rho', @(x,y) 0*x+1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 default_arg('eta', @(x,y) 0*x+1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 dim = 2;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56 C_default = cell(dim,dim,dim,dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 C_default{i,j,k,l} = @(x,y) 0*x ;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 end
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 default_arg('C', C_default);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 assert(isa(g, 'grid.Curvilinear'));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70 if isa(rho, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 rho = grid.evalOn(g, rho);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 if isa(eta, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 eta = grid.evalOn(g, eta);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 C_mat = cell(dim,dim,dim,dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 if isa(C{i,j,k,l}, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 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
85 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 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
87 end
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 obj.C = C_mat;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 elasticObj = scheme.Elastic2dCurvilinearAnisotropic(g, order, rho, C);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95 % Construct a pair of first derivatives
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 K = elasticObj.K;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 K{i,j} = spdiag(K{i,j});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 end
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 J = elasticObj.J;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 Ji = elasticObj.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 D_ref = elasticObj.refObj.D1;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 D1 = cell(dim, 1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 D1Tilde = cell(dim, 1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 D1{i} = 0*D_ref{i};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 D1Tilde{i} = 0*D_ref{i};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 D1{i} = D1{i} + K{i,j}*D_ref{j};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 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
113 end
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 obj.D1 = D1;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 obj.D1Tilde = D1Tilde;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118 eU = elasticObj.E;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 % Storage order for gamma: 11-12-21-22.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 I = speye(g.N(), g.N());
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 eGamma = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 e = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 e{1,1} = [1;0;0;0];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 e{1,2} = [0;1;0;0];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 e{2,1} = [0;0;1;0];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 e{2,2} = [0;0;0;1];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 eGamma{i,j} = kron(I, e{i,j});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 end
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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134 % Store u first, then gamma
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135 mU = dim*g.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 mGamma = dim^2*g.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 Iu = speye(mU, mU);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 Igamma = speye(mGamma, mGamma);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140 Eu = cell2mat({Iu, sparse(mU, mGamma)})';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141 Egamma = cell2mat({sparse(mGamma, mU), Igamma})';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144 eU{i} = Eu*eU{i};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148 eGamma{i,j} = Egamma*eGamma{i,j};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149 end
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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
152 obj.eGamma = eGamma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 obj.eU = eU;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154 obj.Egamma = Egamma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155 obj.Eu = Eu;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
156
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
157 % Build stress operator
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
158 sigma = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 C = obj.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
160 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
161 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
162 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
163 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
164 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
165 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
166 end
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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
171 % Elastic operator
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
172 Delastic = Eu*elasticObj.D*Eu';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174 % Add viscous strains to momentum balance
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 RHOi = spdiag(1./rho);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
176 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
177 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181 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
182 end
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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187 ETA = spdiag(eta);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
188 DstrainRate = 0*Delastic;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
189 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191 DstrainRate = DstrainRate + eGamma{i,j}*(ETA\sigma{i,j});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192 end
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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195 obj.D = Delastic + Dviscous + DstrainRate;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196 obj.Delastic = Delastic;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
197 obj.Dviscous = Dviscous;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198 obj.DstrainRate = DstrainRate;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 obj.sigma = sigma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201 %---- Set remaining object properties ------
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 obj.RHO = elasticObj.RHO;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 obj.ETA = ETA;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 obj.H = elasticObj.H;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206 obj.n_w = elasticObj.n_w;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207 obj.n_e = elasticObj.n_e;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 obj.n_s = elasticObj.n_s;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 obj.n_n = elasticObj.n_n;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 obj.H_w = elasticObj.H_w;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212 obj.H_e = elasticObj.H_e;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
213 obj.H_s = elasticObj.H_s;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214 obj.H_n = elasticObj.H_n;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 obj.e_scalar_w = elasticObj.e_scalar_w;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 obj.e_scalar_e = elasticObj.e_scalar_e;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
218 obj.e_scalar_s = elasticObj.e_scalar_s;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219 obj.e_scalar_n = elasticObj.e_scalar_n;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221 % Misc.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 obj.elasticObj = elasticObj;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 obj.m = elasticObj.m;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 obj.h = elasticObj.h;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226 obj.order = order;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
227 obj.grid = g;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228 obj.dim = dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 end
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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
233 % 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
234 % 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
235 % 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
236 % 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
237 % on the first component. Can also be e.g.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
238 % {'normal', 'd'} or {'tangential', 't'} for conditions on
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
239 % tangential/normal component.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 % 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
241 % 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
242 % neighbour_boundary is a string specifying which boundary to interface to.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
243
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
244 % For displacement bc:
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
245 % bc = {comp, 'd', dComps},
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
246 % where
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
247 % dComps = vector of components with displacement BC. Default: 1:dim.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248 % 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
249 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
250 default_arg('tuning', 1.0);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
251 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
252
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253 component = bc{1};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
254 type = bc{2};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 dim = obj.dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
257 n = obj.getNormal(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
258 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
259 e = obj.getBoundaryOperatorForScalarField('e', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
260
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
261 H = obj.H;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
262 RHO = obj.RHO;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
263 C = obj.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
264 Eu = obj.Eu;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
265 eU = obj.eU;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
266 eGamma = obj.eGamma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
267
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
268 switch type
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
269 case {'F','f','Free','free','traction','Traction','t','T'}
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
271 closure = Eu*closure*Eu';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
272 penalty = Eu*penalty;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
273
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
274 j = component;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
275 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
276 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
277 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
278 closure = closure + eU{j}*( (RHO*H)\(C{i,j,k,l}*e*H_gamma*n{i}*e'*eGamma{k,l}') );
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
279 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
280 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
281 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
282 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
283
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
284 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
285
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
286 function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
287 disp('WARNING: displacementBCNormalTangential is only guaranteed to work for displacement BC on one component and traction on the other');
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
288 u = obj;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
289
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
290 component = bc{1};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
291 type = bc{2};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
292
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
293 switch component
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
294 case 'n'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
295 en = u.getBoundaryOperator('en', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
296 tau_n = u.getBoundaryOperator('tau_n', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
297 N = u.getNormal(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
298 case 't'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
299 en = u.getBoundaryOperator('et', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
300 tau_n = u.getBoundaryOperator('tau_t', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
301 N = u.getTangent(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
302 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
303
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
304 % Operators
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
305 e = u.getBoundaryOperatorForScalarField('e', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
306 h11 = u.getBorrowing(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
307 n = u.getNormal(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
308
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
309 C = u.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
310 Ji = u.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
311 s = spdiag(u.(['s_' boundary]));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
312 m_tot = u.grid.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
313
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
314 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
315 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
316
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
317 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
318 dim = u.dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
319
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
320 % Preallocate
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
321 [~, m_int] = size(H_gamma);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
322 closure = sparse(dim*m_tot, dim*m_tot);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
323 penalty = sparse(dim*m_tot, m_int);
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 % Term 1: The symmetric term
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
326 Z = sparse(m_int, m_int);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
327 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
328 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
329 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
330 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
331 Z = Z + n{i}*N{j}*e'*Ji*C{i,j,k,l}*e*n{k}*N{l};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
332 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
333 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
334 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
335 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
336
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
337 Z = -tuning*dim*1/h11*s*Z;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
338 closure = closure + en*H_gamma*Z*en';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
339 penalty = penalty - en*H_gamma*Z;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
340
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
341 % Term 2: The symmetrizing term
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
342 closure = closure + tau_n*H_gamma*en';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
343 penalty = penalty - tau_n*H_gamma;
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 % Multiply all normal component terms by inverse of density x quadraure
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
346 closure = RHOi*Hi*closure;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
347 penalty = RHOi*Hi*penalty;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
348 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
349
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
350 % type Struct that specifies the interface coupling.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 % Fields:
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352 % -- tuning: penalty strength, defaults to 1.0
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
353 % -- interpolation: type of interpolation, default 'none'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
354 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
355
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
356 defaultType.tuning = 1.0;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
357 defaultType.interpolation = 'none';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
358 defaultType.type = 'standard';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
359 default_struct('type', defaultType);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
360
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
361 switch type.type
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
362 case 'standard'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
363 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
364 case 'frictionalFault'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
365 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
366 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
367
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
368 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
369
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
370 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
371 tuning = type.tuning;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
372
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
373 % u denotes the solution in the own domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
374 % v denotes the solution in the neighbour domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
375
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
376 u = obj;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377 v = neighbour_scheme;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
378
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
379 % Operators, u side
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
380 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
381 en_u = u.getBoundaryOperator('en', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
382 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
383 h11_u = u.getBorrowing(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
384 n_u = u.getNormal(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
385
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
386 C_u = u.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
387 Ji_u = u.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
388 s_u = spdiag(u.(['s_' boundary]));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
389 m_tot_u = u.grid.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
390
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
391 % Operators, v side
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
392 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
393 en_v = v.getBoundaryOperator('en', neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
394 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
395 h11_v = v.getBorrowing(neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
396 n_v = v.getNormal(neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
397
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
398 C_v = v.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
399 Ji_v = v.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
400 s_v = spdiag(v.(['s_' neighbour_boundary]));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
401 m_tot_v = v.grid.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
402
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
403 % Operators that are only required for own domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
404 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
405 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
406
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
407 % Shared operators
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
408 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
409 dim = u.dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
410
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
411 % Preallocate
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
412 [~, m_int] = size(H_gamma);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
413 closure = sparse(dim*m_tot_u, dim*m_tot_u);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
414 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
415
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
416 % Continuity of normal displacement, term 1: The symmetric term
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
417 Z_u = sparse(m_int, m_int);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
418 Z_v = sparse(m_int, m_int);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
419 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
420 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
421 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
422 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
423 Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*n_u{l};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
424 Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*n_v{l};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
425 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
426 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
427 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
428 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
429
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
430 Z = -tuning*dim*( 1/(4*h11_u)*s_u*Z_u + 1/(4*h11_v)*s_v*Z_v );
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
431 closure = closure + en_u*H_gamma*Z*en_u';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
432 penalty = penalty + en_u*H_gamma*Z*en_v';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
433
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
434 % Continuity of normal displacement, term 2: The symmetrizing term
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
435 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
436 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
437
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
438 % Continuity of normal traction
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
439 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
440 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
441
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
442 % Multiply all normal component terms by inverse of density x quadraure
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
443 closure = RHOi*Hi*closure;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
444 penalty = RHOi*Hi*penalty;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
445
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
446 % ---- Tangential tractions are imposed just like traction BC ------
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
447 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
448
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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
451
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
452 % Returns h11 for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
453 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
454 function h11 = getBorrowing(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
455 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 switch boundary
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
458 case {'w','e'}
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
459 h11 = obj.refObj.h11{1};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
460 case {'s', 'n'}
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
461 h11 = obj.refObj.h11{2};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
462 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
463 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
464
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
465 % 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
466 % 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
467 function n = getNormal(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
468 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
469
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
470 n = obj.(['n_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
471 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
472
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
473 % 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
474 % 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
475 function t = getTangent(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
476 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
477
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
478 t = obj.(['tangent_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
479 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
480
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
481 % 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
482 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
483 function o = getBoundaryOperator(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
484 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
485 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
486
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
487 o = obj.([op, '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
488
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
489 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
490
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
491 % 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
492 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
493 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
494 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
495 assertIsMember(op, {'e'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
496
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
497 switch op
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
498
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
499 case 'e'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
500 o = obj.(['e_scalar', '_', 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 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
504
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
505 % 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
506 % Formula: tau_i = T_ij u_j
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
507 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
508 function T = getBoundaryTractionOperator(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
509 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 T = obj.(['T', '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
512 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
513
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
514 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
515 % corresponding to the number of boundary unknowns
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
516 %
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
517 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
518 function H = getBoundaryQuadrature(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
519 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 H = obj.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
522 I_dim = speye(obj.dim, obj.dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
523 H = kron(H, I_dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
524 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
525
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
526 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
527 % corresponding to the number of boundary grid points
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 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
530 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
531 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 H_b = obj.(['H_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
534 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
535
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
536 function N = size(obj)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
537 N = (obj.dim + obj.dim^2)*prod(obj.m);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
538 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
539 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
540 end