Mercurial > repos > public > sbplib

annotate +scheme/ViscoElastic2d.m @ 1310:eb015fe9605b feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
author Martin Almquist <malmquist@stanford.edu>
date Mon, 20 Jul 2020 22:06:33 -0700
parents f59b849df30f
children b2e84fe336d8
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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
275 switch type
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
276 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
277
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
278 % Get elastic closure and penalty
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
279 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
280 closure = Eu*closure*Eu';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
281 penalty = Eu*penalty;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
282
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
283 switch component
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
284 case 't'
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
285 dir = obj.getTangent(boundary);
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
286 case 'n'
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
287 dir = obj.getNormal(boundary);
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
288 case 1
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
289 dir = {1, 0};
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
290 case 2
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
291 dir = {0, 1};
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
292 end
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 % Get elastic closure and penalty
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
310 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning);
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
311 closure = Eu*closure*Eu';
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
312 penalty = Eu*penalty;
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
313
1310
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
314 switch component
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
315 case 't'
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
316 dir = obj.getTangent(boundary);
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
317 case 'n'
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
318 dir = obj.getNormal(boundary);
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
319 case 1
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
320 dir = {1, 0};
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
321 case 2
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
322 dir = {0, 1};
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
323 end
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
324
1309
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
325 % Add penalty to strain rate eq
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
326 for i = 1:dim
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
327 for j = 1:dim
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
328 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
329 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
330 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
331 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
332 end
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
333 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
334 end
1309
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
335 end
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
336 end
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
337 end
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
338
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
339 end
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 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
342
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
343 % type Struct that specifies the interface coupling.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
344 % Fields:
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
345 % -- tuning: penalty strength, defaults to 1.0
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
346 % -- interpolation: type of interpolation, default 'none'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
347 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
348
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
349 defaultType.tuning = 1.0;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
350 defaultType.interpolation = 'none';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 defaultType.type = 'standard';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352 default_struct('type', defaultType);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
353
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
354 switch type.type
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
355 case 'standard'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
356 [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
357 case 'frictionalFault'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
358 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
359 end
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 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
362
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
363 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
364 tuning = type.tuning;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
365
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
366 % u denotes the solution in the own domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
367 % v denotes the solution in the neighbour domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
368
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
369 u = obj;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
370 v = neighbour_scheme;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
371
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
372 % Operators, u side
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
373 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
374 en_u = u.getBoundaryOperator('en', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
375 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
376 h11_u = u.getBorrowing(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377 n_u = u.getNormal(boundary);
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 C_u = u.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
380 Ji_u = u.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
381 s_u = spdiag(u.(['s_' boundary]));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
382 m_tot_u = u.grid.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
383
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
384 % Operators, v side
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
385 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
386 en_v = v.getBoundaryOperator('en', neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
387 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
388 h11_v = v.getBorrowing(neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
389 n_v = v.getNormal(neighbour_boundary);
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 C_v = v.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
392 Ji_v = v.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
393 s_v = spdiag(v.(['s_' neighbour_boundary]));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
394 m_tot_v = v.grid.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
395
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
396 % Operators that are only required for own domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
397 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
398 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
399
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
400 % Shared operators
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
401 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
402 dim = u.dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
403
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
404 % Preallocate
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
405 [~, m_int] = size(H_gamma);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
406 closure = sparse(dim*m_tot_u, dim*m_tot_u);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
407 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
408
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
409 % Continuity of normal displacement, term 1: The symmetric term
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
410 Z_u = sparse(m_int, m_int);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
411 Z_v = sparse(m_int, m_int);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
412 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
413 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
414 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
415 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
416 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
417 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
418 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
419 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
420 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
421 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
422
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
423 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
424 closure = closure + en_u*H_gamma*Z*en_u';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
425 penalty = penalty + en_u*H_gamma*Z*en_v';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
426
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
427 % Continuity of normal displacement, term 2: The symmetrizing term
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
428 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
429 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
430
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
431 % Continuity of normal traction
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
432 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
433 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
434
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
435 % Multiply all normal component terms by inverse of density x quadraure
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
436 closure = RHOi*Hi*closure;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
437 penalty = RHOi*Hi*penalty;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
438
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
439 % ---- Tangential tractions are imposed just like traction BC ------
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
440 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
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 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
443
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
444
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
445 % Returns h11 for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
446 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
447 function h11 = getBorrowing(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
448 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
449
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
450 switch boundary
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
451 case {'w','e'}
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
452 h11 = obj.refObj.h11{1};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
453 case {'s', 'n'}
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
454 h11 = obj.refObj.h11{2};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
455 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
456 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
457
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
458 % 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
459 % 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
460 function n = getNormal(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
461 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
462
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
463 n = obj.(['n_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
464 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
465
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
466 % 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
467 % 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
468 function t = getTangent(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
469 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
470
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
471 t = obj.(['tangent_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
472 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
473
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
474 % 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
475 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
476 function o = getBoundaryOperator(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
477 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
478 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
479
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
480 o = obj.([op, '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
481
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
482 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
483
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
484 % 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
485 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
486 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
487 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
488 assertIsMember(op, {'e'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
489
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
490 switch op
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
491
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
492 case 'e'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
493 o = obj.(['e_scalar', '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
494 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
495
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
496 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
497
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
498 % 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
499 % Formula: tau_i = T_ij u_j
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
500 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
501 function T = getBoundaryTractionOperator(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
502 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
503
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
504 T = obj.(['T', '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
505 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
506
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
507 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
508 % corresponding to the number of boundary unknowns
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
509 %
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
510 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
511 function H = getBoundaryQuadrature(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
512 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 H = obj.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
515 I_dim = speye(obj.dim, obj.dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
516 H = kron(H, I_dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
517 end
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 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
520 % corresponding to the number of boundary grid points
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
521 %
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
522 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
523 function H_b = getBoundaryQuadratureForScalarField(obj, 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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
526 H_b = obj.(['H_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
527 end
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 function N = size(obj)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
530 N = (obj.dim + obj.dim^2)*prod(obj.m);
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 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
533 end