Mercurial > repos > public > sbplib

annotate +scheme/ViscoElastic2d.m @ 1334:df8c71b80c33 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use the logic grid associated with a CurvilinearGrid instead of creating a new one
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sat, 07 May 2022 10:40:47 +0200
parents cd0bfcaef0ba
children
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
1314
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
45 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
46 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
47 tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
48 tau_t_w, tau_t_e, tau_t_s, tau_t_n % Return scalar field
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
49
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 elasticObj
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 methods
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56 % The coefficients can either be function handles or grid functions
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 function obj = ViscoElastic2d(g, order, rho, C, eta)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 default_arg('rho', @(x,y) 0*x+1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 default_arg('eta', @(x,y) 0*x+1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 dim = 2;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 C_default = cell(dim,dim,dim,dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 C_default{i,j,k,l} = @(x,y) 0*x ;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 default_arg('C', C_default);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 assert(isa(g, 'grid.Curvilinear'));
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 if isa(rho, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 rho = grid.evalOn(g, rho);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 if isa(eta, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 eta = grid.evalOn(g, eta);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 C_mat = cell(dim,dim,dim,dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 if isa(C{i,j,k,l}, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 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
91 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 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
93 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 obj.C = C_mat;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 elasticObj = scheme.Elastic2dCurvilinearAnisotropic(g, order, rho, C);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101 % Construct a pair of first derivatives
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 K = elasticObj.K;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 K{i,j} = spdiag(K{i,j});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 J = elasticObj.J;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 Ji = elasticObj.Ji;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 D_ref = elasticObj.refObj.D1;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 D1 = cell(dim, 1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 D1Tilde = cell(dim, 1);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 D1{i} = 0*D_ref{i};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 D1Tilde{i} = 0*D_ref{i};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 D1{i} = D1{i} + K{i,j}*D_ref{j};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118 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
119 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 obj.D1 = D1;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 obj.D1Tilde = D1Tilde;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 eU = elasticObj.E;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 % Storage order for gamma: 11-12-21-22.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 I = speye(g.N(), g.N());
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128 eGamma = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129 e = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 e{1,1} = [1;0;0;0];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 e{1,2} = [0;1;0;0];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132 e{2,1} = [0;0;1;0];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 e{2,2} = [0;0;0;1];
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 eGamma{i,j} = kron(I, e{i,j});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140 % Store u first, then gamma
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141 mU = dim*g.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142 mGamma = dim^2*g.N();
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143 Iu = speye(mU, mU);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144 Igamma = speye(mGamma, mGamma);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 Eu = cell2mat({Iu, sparse(mU, mGamma)})';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 Egamma = cell2mat({sparse(mGamma, mU), Igamma})';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
150 eU{i} = Eu*eU{i};
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 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154 eGamma{i,j} = Egamma*eGamma{i,j};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
156 end
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 obj.eGamma = eGamma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 obj.eU = eU;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
160 obj.Egamma = Egamma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
161 obj.Eu = Eu;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
162
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
163 % Build stress operator
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
164 sigma = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
165 C = obj.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
166 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
167 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168 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
169 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
170 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
171 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
172 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
176
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
177 % Elastic operator
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178 Delastic = Eu*elasticObj.D*Eu';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 % Add viscous strains to momentum balance
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181 RHOi = spdiag(1./rho);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
182 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
183 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
185 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187 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
188 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
189 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
193 ETA = spdiag(eta);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
194 DstrainRate = 0*Delastic;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
197 DstrainRate = DstrainRate + eGamma{i,j}*(ETA\sigma{i,j});
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201 obj.D = Delastic + Dviscous + DstrainRate;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 obj.Delastic = Delastic;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 obj.Dviscous = Dviscous;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 obj.DstrainRate = DstrainRate;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 obj.sigma = sigma;
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 %---- Set remaining object properties ------
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 obj.RHO = elasticObj.RHO;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 obj.ETA = ETA;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210 obj.H = elasticObj.H;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212 obj.n_w = elasticObj.n_w;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
213 obj.n_e = elasticObj.n_e;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214 obj.n_s = elasticObj.n_s;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 obj.n_n = elasticObj.n_n;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
217 obj.tangent_w = elasticObj.tangent_w;
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
218 obj.tangent_e = elasticObj.tangent_e;
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
219 obj.tangent_s = elasticObj.tangent_s;
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
220 obj.tangent_n = elasticObj.tangent_n;
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
221
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 obj.H_w = elasticObj.H_w;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 obj.H_e = elasticObj.H_e;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 obj.H_s = elasticObj.H_s;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 obj.H_n = elasticObj.H_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 obj.e_scalar_w = elasticObj.e_scalar_w;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228 obj.e_scalar_e = elasticObj.e_scalar_e;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 obj.e_scalar_s = elasticObj.e_scalar_s;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 obj.e_scalar_n = elasticObj.e_scalar_n;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231
1314
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
232 % -- Create new traction operators including viscous strain contribution --
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
233 tau1 = struct;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
234 tau2 = struct;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
235 tau_n = struct;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
236 tau_t = struct;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
237 boundaries = {'w', 'e', 's', 'n'};
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
238 for bNumber = 1:numel(boundaries)
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
239 b = boundaries{bNumber};
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
240
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
241 n = elasticObj.getNormal(b);
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
242 t = elasticObj.getTangent(b);
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
243 e = elasticObj.getBoundaryOperatorForScalarField('e', b);
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
244 tau1.(b) = (elasticObj.getBoundaryOperator('tau1', b)'*Eu')';
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
245 tau2.(b) = (elasticObj.getBoundaryOperator('tau2', b)'*Eu')';
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
246
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
247 % Add viscous contributions
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
248 for i = 1:dim
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
249 for k = 1:dim
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
250 for l = 1:dim
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
251 tau1.(b) = tau1.(b) - (n{i}*e'*C{i,1,k,l}*eGamma{k,l}')';
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
252 tau2.(b) = tau2.(b) - (n{i}*e'*C{i,2,k,l}*eGamma{k,l}')';
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
253 end
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
254 end
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
255 end
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
256
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
257 % Compute normal and tangential components
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
258 tau_n.(b) = tau1.(b)*n{1} + tau2.(b)*n{2};
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
259 tau_t.(b) = tau1.(b)*t{1} + tau2.(b)*t{2};
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
260 end
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
261
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
262 obj.tau1_w = tau1.w;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
263 obj.tau1_e = tau1.e;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
264 obj.tau1_s = tau1.s;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
265 obj.tau1_n = tau1.n;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
266
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
267 obj.tau2_w = tau2.w;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
268 obj.tau2_e = tau2.e;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
269 obj.tau2_s = tau2.s;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
270 obj.tau2_n = tau2.n;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
271
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
272 obj.tau_n_w = tau_n.w;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
273 obj.tau_n_e = tau_n.e;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
274 obj.tau_n_s = tau_n.s;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
275 obj.tau_n_n = tau_n.n;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
276
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
277 obj.tau_t_w = tau_t.w;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
278 obj.tau_t_e = tau_t.e;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
279 obj.tau_t_s = tau_t.s;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
280 obj.tau_t_n = tau_t.n;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
281 %----------------------------------------
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
282
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
283 % Misc.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
284 obj.elasticObj = elasticObj;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
285 obj.m = elasticObj.m;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
286 obj.h = elasticObj.h;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
287
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
288 obj.order = order;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
289 obj.grid = g;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
290 obj.dim = dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
291
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
292 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
293
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
294
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
295 % 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
296 % 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
297 % 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
298 % 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
299 % on the first component. Can also be e.g.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
300 % {'normal', 'd'} or {'tangential', 't'} for conditions on
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
301 % tangential/normal component.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
302 % 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
303 % 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
304 % neighbour_boundary is a string specifying which boundary to interface to.
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 % For displacement bc:
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
307 % bc = {comp, 'd', dComps},
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
308 % where
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
309 % dComps = vector of components with displacement BC. Default: 1:dim.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
310 % 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
311 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
312 default_arg('tuning', 1.0);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
313 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
314
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
315 component = bc{1};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
316 type = bc{2};
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
317 dim = obj.dim;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
318
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
319 n = obj.getNormal(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
320 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
321 e = obj.getBoundaryOperatorForScalarField('e', boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
322
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
323 H = obj.H;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
324 RHO = obj.RHO;
1309
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
325 ETA = obj.ETA;
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
326 C = obj.C;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
327 Eu = obj.Eu;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
328 eU = obj.eU;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
329 eGamma = obj.eGamma;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
330
1311
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
331 % Get elastic closure and penalty
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
332 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning);
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
333 closure = Eu*closure*Eu';
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
334 penalty = Eu*penalty;
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
335
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
336 switch component
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
337 case 't'
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
338 dir = obj.getTangent(boundary);
1314
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
339 tau = obj.getBoundaryOperator('tau_t', boundary);
1311
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
340 case 'n'
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
341 dir = obj.getNormal(boundary);
1314
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
342 tau = obj.getBoundaryOperator('tau_n', boundary);
1311
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
343 case 1
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
344 dir = {1, 0};
1314
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
345 tau = obj.getBoundaryOperator('tau1', boundary);
1311
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
346 case 2
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
347 dir = {0, 1};
1314
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
348 tau = obj.getBoundaryOperator('tau2', boundary);
1311
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
349 end
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement
Martin Almquist <malmquist@stanford.edu>
parents: 1310
diff changeset
350
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 switch type
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352 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
353
1314
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
354 % Set elastic closure to zero
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
355 closure = 0*closure;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
356
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
357 for m = 1:dim
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
358 closure = closure - eU{m}*( (RHO*H)\(e*dir{m}*H_gamma*tau') );
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
359 end
1308
5016f3f3badb Add viscoElastic traction bc for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1307
diff changeset
360
1309
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
361 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
362
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
363 % Add penalty to strain rate eq
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
364 for i = 1:dim
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
365 for j = 1:dim
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
366 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
367 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
368 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
369 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
370 end
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential
Martin Almquist <malmquist@stanford.edu>
parents: 1309
diff changeset
371 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
372 end
1309
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
373 end
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
374 end
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
375 end
f59b849df30f Add displacement BC to viscoelastic
Martin Almquist <malmquist@stanford.edu>
parents: 1308
diff changeset
376
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377 end
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 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
380
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
381 % type Struct that specifies the interface coupling.
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
382 % Fields:
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
383 % -- tuning: penalty strength, defaults to 1.0
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
384 % -- interpolation: type of interpolation, default 'none'
1320
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
385 function [closure, penalty, forcingPenalties] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
386
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
387 defaultType.tuning = 1.0;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
388 defaultType.interpolation = 'none';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
389 defaultType.type = 'standard';
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
390 default_struct('type', defaultType);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
391
1320
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
392 forcingPenalties = [];
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
393
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
394 switch type.type
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
395 case 'standard'
1312
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
396 [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
397 case 'frictionalFault'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
398 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
1318
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
399 case 'normalTangential'
1320
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
400 [closure, penalty, forcingPenalties] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
401 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
402
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
403 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
404
1312
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
405 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
406
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
407 % u denotes the solution in the own domain
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
408 % v denotes the solution in the neighbour domain
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
409 u = obj;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
410 v = neighbour_scheme;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
411
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
412 dim = obj.dim;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
413
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
414 n = u.getNormal(boundary);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
415 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
416 e = u.getBoundaryOperatorForScalarField('e', boundary);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
417
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
418 ev = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
419
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
420 H = u.H;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
421 RHO = u.RHO;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
422 ETA = u.ETA;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
423 C = u.C;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
424 Eu = u.Eu;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
425 eU = u.eU;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
426 eGamma = u.eGamma;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
427
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
428 CV = v.C;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
429 Ev = v.Eu;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
430 eV = v.eU;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
431 eGammaV = v.eGamma;
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
432 nV = v.getNormal(neighbour_boundary);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
433
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
434
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
435 % Get elastic closure and penalty
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
436 [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
437 closure = Eu*closure*Eu';
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
438 penalty = Eu*penalty*Ev';
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
439
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
440 % Add viscous part of traction coupling
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
441 for i = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
442 for j = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
443 for k = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
444 for l = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
445 closure = closure + 1/2*eU{j}*( (RHO*H)\(C{i,j,k,l}*e*H_gamma*n{i}*e'*eGamma{k,l}') );
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
446 penalty = penalty + 1/2*eU{j}*( (RHO*H)\(e*H_gamma*nV{i}*(ev'*CV{i,j,k,l}*ev)*ev'*eGammaV{k,l}') );
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
447 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
448 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
449 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
450 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
451
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
452 % Add penalty to strain rate eq
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
453 for i = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
454 for j = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
455 for k = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
456 for l = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
457 closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}*e'*eU{l}') );
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
458 penalty = penalty + 1/2*eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}*ev'*eV{l}') );
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
459 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
460 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
461 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
462 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
463
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
464
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
465 end
6d64b57caf46 Add viscoelastic regular interfaces
Martin Almquist <malmquist@stanford.edu>
parents: 1311
diff changeset
466
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
467 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
468 tuning = type.tuning;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
469
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
470 % u denotes the solution in the own domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
471 % v denotes the solution in the neighbour domain
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
472 u = obj;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
473 v = neighbour_scheme;
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
474
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
475 dim = obj.dim;
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
476
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
477 n = u.getNormal(boundary);
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
478 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
479 e = u.getBoundaryOperatorForScalarField('e', boundary);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
480
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
481 ev = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
482
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
483 H = u.H;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
484 RHO = u.RHO;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
485 ETA = u.ETA;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
486 C = u.C;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
487 Eu = u.Eu;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
488 eU = u.eU;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
489 eGamma = u.eGamma;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
490 Egamma = u.Egamma;
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
491
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
492 CV = v.C;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
493 Ev = v.Eu;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
494 eV = v.eU;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
495 eGammaV = v.eGamma;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
496 nV = v.getNormal(neighbour_boundary);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
497
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
498 % Reduce stiffness tensors to boundary size
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
499 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
500 for j = 1:dim
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
501 for k = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
502 for l = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
503 C{i,j,k,l} = e'*C{i,j,k,l}*e;
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
504 CV{i,j,k,l} = ev'*CV{i,j,k,l}*ev;
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
505 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
506 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
507 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
508 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
509
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
510 % Get elastic closure and penalty
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
511 [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
512 closure = Eu*closure*Eu';
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
513 penalty = Eu*penalty*Ev';
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
514
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
515 % ---- Tangential tractions are imposed just like traction BC ------
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
516 % We only need the viscous part
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
517 closure_tangential = obj.boundary_condition(boundary, {'t', 't'});
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
518 closure = closure + closure_tangential*Egamma*Egamma';
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
519
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
520
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
521 % ------ Coupling of normal component -----------
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
522 % Add viscous part of traction coupling
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
523 for i = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
524 for j = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
525 for k = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
526 for l = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
527 for m = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
528 closure = closure + 1/2*eU{m}*( (RHO*H)\(e*n{m}*H_gamma*n{j}*n{i}*C{i,j,k,l}*e'*eGamma{k,l}') );
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
529 penalty = penalty - 1/2*eU{m}*( (RHO*H)\(e*n{m}*H_gamma*nV{j}*nV{i}*CV{i,j,k,l}*ev'*eGammaV{k,l}') );
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
530 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
531 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
532 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
533 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
534 end
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
535
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
536 % Add penalty to strain rate eq
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
537 for i = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
538 for j = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
539 for k = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
540 for l = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
541 for m = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
542 closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma*n{m}*e'*eU{m}') );
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
543 penalty = penalty - 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma*nV{m}*ev'*eV{m}') );
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
544 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
545 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
546 end
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
547 end
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
548 end
1313
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
549 %-------------------------------------------------
a4c8bf2371f3 Add viscoelastic fault interface
Martin Almquist <malmquist@stanford.edu>
parents: 1312
diff changeset
550
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
551 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
552
1320
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
553 function [closure, penalty, forcingPenalties] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type)
1318
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
554 tuning = type.tuning;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
555
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
556 % u denotes the solution in the own domain
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
557 % v denotes the solution in the neighbour domain
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
558 u = obj;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
559 v = neighbour_scheme;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
560
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
561 dim = obj.dim;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
562
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
563 n = u.getNormal(boundary);
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
564 t = u.getTangent(boundary);
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
565 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
566 e = u.getBoundaryOperatorForScalarField('e', boundary);
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
567
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
568 ev = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
569
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
570 H = u.H;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
571 RHO = u.RHO;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
572 ETA = u.ETA;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
573 C = u.C;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
574 Eu = u.Eu;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
575 eU = u.eU;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
576 eGamma = u.eGamma;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
577 Egamma = u.Egamma;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
578
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
579 CV = v.C;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
580 Ev = v.Eu;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
581 eV = v.eU;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
582 eGammaV = v.eGamma;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
583 nV = v.getNormal(neighbour_boundary);
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
584 tV = v.getTangent(neighbour_boundary);
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
585
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
586 % Reduce stiffness tensors to boundary size
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
587 for i = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
588 for j = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
589 for k = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
590 for l = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
591 C{i,j,k,l} = e'*C{i,j,k,l}*e;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
592 CV{i,j,k,l} = ev'*CV{i,j,k,l}*ev;
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
593 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
594 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
595 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
596 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
597
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
598 % Get elastic closure and penalty
1320
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
599 [closure, penalty, forcingPenalties] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
1318
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
600 closure = Eu*closure*Eu';
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
601 penalty = Eu*penalty*Ev';
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
602
1320
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
603 for i = 1:numel(forcingPenalties)
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
604 forcingPenalties{i} = Eu*forcingPenalties{i};
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
605 end
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
606 forcing_u_n = forcingPenalties{1};
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
607 forcing_u_t = forcingPenalties{2};
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
608
1318
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
609 % ------ Traction coupling, viscous part -----------
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
610 for i = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
611 for j = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
612 for k = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
613 for l = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
614 for m = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
615 % Normal component
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
616 closure = closure + 1/2*eU{m}*( (RHO*H)\(e*n{m}*H_gamma*n{j}*n{i}*C{i,j,k,l}*e'*eGamma{k,l}') );
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
617 penalty = penalty - 1/2*eU{m}*( (RHO*H)\(e*n{m}*H_gamma*nV{j}*nV{i}*CV{i,j,k,l}*ev'*eGammaV{k,l}') );
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
618
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
619 % Tangential component
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
620 closure = closure + 1/2*eU{m}*( (RHO*H)\(e*t{m}*H_gamma*t{j}*n{i}*C{i,j,k,l}*e'*eGamma{k,l}') );
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
621 penalty = penalty - 1/2*eU{m}*( (RHO*H)\(e*t{m}*H_gamma*tV{j}*nV{i}*CV{i,j,k,l}*ev'*eGammaV{k,l}') );
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
622 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
623 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
624 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
625 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
626 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
627 %-------------------------------------------------
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
628
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
629 % --- Displacement coupling ----------------------
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
630 % Add penalty to strain rate eq
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
631 for i = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
632 for j = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
633 for k = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
634 for l = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
635 for m = 1:dim
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
636 % Normal component
1320
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
637 closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma*n{m}*e'*eU{m}') );
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
638 penalty = penalty - 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma*nV{m}*ev'*eV{m}') );
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
639
1318
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
640
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
641 % Tangential component
1320
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
642 closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(e*t{l}*n{k}*C{i,j,k,l}*H_gamma*t{m}*e'*eU{m}') );
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
643 penalty = penalty - 1/2*eGamma{i,j}*( (H*ETA)\(e*t{l}*n{k}*C{i,j,k,l}*H_gamma*tV{m}*ev'*eV{m}') );
1318
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
644 end
1320
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
645 forcing_u_n = forcing_u_n + 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma) );
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
646 forcing_u_t = forcing_u_t + 1/2*eGamma{i,j}*( (H*ETA)\(e*t{l}*n{k}*C{i,j,k,l}*H_gamma) );
1318
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
647 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
648 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
649 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
650 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
651 %-------------------------------------------------
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
652
1320
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
653 forcingPenalties{1} = forcing_u_n;
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
654 forcingPenalties{2} = forcing_u_t;
cd0bfcaef0ba Add some interface forcing penalties in ViscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1318
diff changeset
655
1318
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
656 end
6650b111742d Implement normal-tangential interface treatment in viscoElastic2d
Martin Almquist <malmquist@stanford.edu>
parents: 1314
diff changeset
657
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
658 % 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
659 % 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
660 function n = getNormal(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
661 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
662
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
663 n = obj.(['n_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
664 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
665
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
666 % 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
667 % 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
668 function t = getTangent(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
669 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
670
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
671 t = obj.(['tangent_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
672 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
673
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
674 % 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
675 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
676 function o = getBoundaryOperator(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
677 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
678 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
679
1314
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
Martin Almquist <malmquist@stanford.edu>
parents: 1313
diff changeset
680 o = obj.([op, '_', boundary]);
1307
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
681
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
682 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
683
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
684 % 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
685 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
686 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
687 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
688 assertIsMember(op, {'e'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
689
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
690 switch op
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
691
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
692 case 'e'
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
693 o = obj.(['e_scalar', '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
694 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
695
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
696 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
697
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
698 % 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
699 % Formula: tau_i = T_ij u_j
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
700 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
701 function T = getBoundaryTractionOperator(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
702 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
703
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
704 T = obj.(['T', '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
705 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
706
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
707 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
708 % corresponding to the number of boundary unknowns
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
709 %
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
710 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
711 function H = getBoundaryQuadrature(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
712 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
713
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
714 H = obj.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
715 I_dim = speye(obj.dim, obj.dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
716 H = kron(H, I_dim);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
717 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
718
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
719 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
720 % corresponding to the number of boundary grid points
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
721 %
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
722 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
723 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
724 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
725
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
726 H_b = obj.(['H_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
727 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
728
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
729 function N = size(obj)
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
730 N = (obj.dim + obj.dim^2)*prod(obj.m);
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
731 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
732 end
fcca6ad8b102 Add diffOp for viscoElastic
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
733 end