Mercurial > repos > public > sbplib

annotate +scheme/ViscoElastic2d.m @ 1308:5016f3f3badb feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add viscoElastic traction bc for normal-tangential
author Martin Almquist <malmquist@stanford.edu>
date Sun, 19 Jul 2020 21:24:48 -0700
parents fcca6ad8b102
children f59b849df30f
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;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
269 C = obj.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270 Eu = obj.Eu;
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 eGamma = obj.eGamma;
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 switch type
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
275 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
276
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
277 % Get elastic closure and penalty
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
278 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
279 closure = Eu*closure*Eu';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
280 penalty = Eu*penalty;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
281
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
282 switch component
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
283 case 't'
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
284 dir = obj.getTangent(boundary);
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
285 case 'n'
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
286 dir = obj.getNormal(boundary);
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
287 case 1
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
288 dir = {1, 0};
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
289 case 2
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
290 dir = {0, 1};
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
291 end
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
292
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
293 % Add viscous part of closure
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
294 for j = 1:dim
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
295 for i = 1:dim
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
296 for k = 1:dim
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
297 for l = 1:dim
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
298 closure = closure + eU{j}*( (RHO*H)\(C{i,j,k,l}*e*dir{j}^2*H_gamma*n{i}*e'*eGamma{k,l}') );
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
299 end
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
300 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
301 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
302 end
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
303
1307
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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
306 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
308 function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
309 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
310 u = obj;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
311
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
312 component = bc{1};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
313 type = bc{2};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
314
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
315 switch component
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
316 case 'n'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
317 en = u.getBoundaryOperator('en', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
318 tau_n = u.getBoundaryOperator('tau_n', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
319 N = u.getNormal(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
320 case 't'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
321 en = u.getBoundaryOperator('et', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
322 tau_n = u.getBoundaryOperator('tau_t', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
323 N = u.getTangent(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
324 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
325
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
326 % Operators
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
327 e = u.getBoundaryOperatorForScalarField('e', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
328 h11 = u.getBorrowing(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
329 n = u.getNormal(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
330
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
331 C = u.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
332 Ji = u.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
333 s = spdiag(u.(['s_' boundary]));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
334 m_tot = u.grid.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
335
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
336 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
337 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
338
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
339 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
340 dim = u.dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
341
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
342 % Preallocate
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
343 [~, m_int] = size(H_gamma);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
344 closure = sparse(dim*m_tot, dim*m_tot);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
345 penalty = sparse(dim*m_tot, m_int);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
346
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
347 % Term 1: The symmetric term
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
348 Z = sparse(m_int, m_int);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
349 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
350 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
353 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
354 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
355 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
356 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
357 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
358
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
359 Z = -tuning*dim*1/h11*s*Z;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
360 closure = closure + en*H_gamma*Z*en';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
361 penalty = penalty - en*H_gamma*Z;
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 % Term 2: The symmetrizing term
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
364 closure = closure + tau_n*H_gamma*en';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
365 penalty = penalty - tau_n*H_gamma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
366
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
367 % Multiply all normal component terms by inverse of density x quadraure
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
368 closure = RHOi*Hi*closure;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
369 penalty = RHOi*Hi*penalty;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
370 end
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 % type Struct that specifies the interface coupling.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
373 % Fields:
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
374 % -- tuning: penalty strength, defaults to 1.0
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
375 % -- interpolation: type of interpolation, default 'none'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
376 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
378 defaultType.tuning = 1.0;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
379 defaultType.interpolation = 'none';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
380 defaultType.type = 'standard';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
381 default_struct('type', defaultType);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
382
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
383 switch type.type
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
384 case 'standard'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
385 [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
386 case 'frictionalFault'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
387 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
388 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
389
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
390 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
391
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
392 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
393 tuning = type.tuning;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
394
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
395 % u denotes the solution in the own domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
396 % v denotes the solution in the neighbour domain
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 u = obj;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
399 v = neighbour_scheme;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
400
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
401 % Operators, u side
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
402 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
403 en_u = u.getBoundaryOperator('en', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
404 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
405 h11_u = u.getBorrowing(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
406 n_u = u.getNormal(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
407
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
408 C_u = u.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
409 Ji_u = u.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
410 s_u = spdiag(u.(['s_' boundary]));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
411 m_tot_u = u.grid.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
412
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
413 % Operators, v side
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
414 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
415 en_v = v.getBoundaryOperator('en', neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
416 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
417 h11_v = v.getBorrowing(neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
418 n_v = v.getNormal(neighbour_boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
419
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
420 C_v = v.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
421 Ji_v = v.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
422 s_v = spdiag(v.(['s_' neighbour_boundary]));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
423 m_tot_v = v.grid.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
424
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
425 % Operators that are only required for own domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
426 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
427 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
428
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
429 % Shared operators
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
430 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
431 dim = u.dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
432
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
433 % Preallocate
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
434 [~, m_int] = size(H_gamma);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
435 closure = sparse(dim*m_tot_u, dim*m_tot_u);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
436 penalty = sparse(dim*m_tot_u, dim*m_tot_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 displacement, term 1: The symmetric term
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
439 Z_u = sparse(m_int, m_int);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
440 Z_v = sparse(m_int, m_int);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
441 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
442 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
443 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
444 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
445 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
446 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
447 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
448 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
449 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
450 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
451
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
452 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
453 closure = closure + en_u*H_gamma*Z*en_u';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
454 penalty = penalty + en_u*H_gamma*Z*en_v';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
455
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
456 % Continuity of normal displacement, term 2: The symmetrizing term
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
457 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
458 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
459
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
460 % Continuity of normal traction
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
461 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
462 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
463
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
464 % Multiply all normal component terms by inverse of density x quadraure
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
465 closure = RHOi*Hi*closure;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
466 penalty = RHOi*Hi*penalty;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
467
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
468 % ---- Tangential tractions are imposed just like traction BC ------
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
469 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
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 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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
474 % Returns h11 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 h11 = getBorrowing(obj, 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
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
479 switch boundary
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
480 case {'w','e'}
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
481 h11 = obj.refObj.h11{1};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
482 case {'s', 'n'}
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
483 h11 = obj.refObj.h11{2};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
484 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
485 end
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 % 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
488 % 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
489 function n = getNormal(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
490 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 n = obj.(['n_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
493 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
494
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
495 % Returns the unit tangent vector for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
496 % 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
497 function t = getTangent(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
498 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
499
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
500 t = obj.(['tangent_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
501 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
502
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
503 % Returns the boundary operator op for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
504 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
505 function o = getBoundaryOperator(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
506 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
507 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
508
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
509 o = obj.([op, '_', boundary]);
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 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
512
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
513 % 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
514 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
515 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
516 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
517 assertIsMember(op, {'e'})
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 switch op
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 case 'e'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
522 o = obj.(['e_scalar', '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
523 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
524
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
525 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
526
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
527 % 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
528 % Formula: tau_i = T_ij u_j
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
529 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
530 function T = getBoundaryTractionOperator(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 T = obj.(['T', '_', 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 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
537 % corresponding to the number of boundary unknowns
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
538 %
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
539 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
540 function H = getBoundaryQuadrature(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
541 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
542
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
543 H = obj.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
544 I_dim = speye(obj.dim, obj.dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
545 H = kron(H, I_dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
546 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
547
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
548 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
549 % corresponding to the number of boundary grid points
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
550 %
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
551 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
552 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
553 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
554
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
555 H_b = obj.(['H_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
556 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
557
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
558 function N = size(obj)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
559 N = (obj.dim + obj.dim^2)*prod(obj.m);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
560 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
561 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
562 end