Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1306:633757e582e5 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add transformation gradient as property in CurvilinearAnisotropic
author Martin Almquist <malmquist@stanford.edu>
date Sun, 19 Jul 2020 20:29:48 -0700
parents a0d615bde7f8
children 6c308da9dcbc
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 classdef Elastic2dCurvilinearAnisotropic < scheme.Scheme
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 % Discretizes the elastic wave equation:
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5 % in curvilinear coordinates.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 % opSet should be cell array of opSets, one per dimension. This
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 % is useful if we have periodic BC in one direction.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 % Assumes fully compatible operators.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 properties
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 m % Number of points in each direction, possibly a vector
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12 h % Grid spacing
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14 grid
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 order % Order of accuracy for the approximation
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 % Diagonal matrices for variable coefficients
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 J, Ji
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21 RHO % Density
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 C % Elastic stiffness tensor
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24 D % Total operator
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25
1306
633757e582e5 Add transformation gradient as property in CurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1302
diff changeset
26 K % Transformation gradient
1211
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
27 Dx, Dy % Physical derivatives
1214
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
28 sigma % Cell matrix of physical stress operators
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
29 n_w, n_e, n_s, n_n % Physical normals
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
30 tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents
1211
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
31
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32 % Boundary operators in cell format, used for BC
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 T_w, T_e, T_s, T_n
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35 % Traction operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36 tau_w, tau_e, tau_s, tau_n % Return vector field
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
39 tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
40 tau_t_w, tau_t_e, tau_t_s, tau_t_n % Return scalar field
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 % Inner products
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 H
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45 % Boundary inner products (for scalar field)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 H_w, H_e, H_s, H_n
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 % Surface Jacobian vectors
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 s_w, s_e, s_s, s_n
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 % Boundary restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
56 en_w, en_e, en_s, en_n % Act on vector field, return normal component
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
57 et_w, et_e, et_s, et_n % Act on vector field, return tangential component
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 % E{i}^T picks out component i
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 E
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 % Elastic2dVariableAnisotropic object for reference domain
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 refObj
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 methods
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 % The coefficients can either be function handles or grid functions
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
1302
a0d615bde7f8 Add the hollow option to the anisotropic diffops
Martin Almquist <malmquist@stanford.edu>
parents: 1295
diff changeset
70 function obj = Elastic2dCurvilinearAnisotropic(g, order, rho, C, opSet, optFlag, hollow)
a0d615bde7f8 Add the hollow option to the anisotropic diffops
Martin Almquist <malmquist@stanford.edu>
parents: 1295
diff changeset
71 default_arg('hollow', false);
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 default_arg('rho', @(x,y) 0*x+1);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 default_arg('optFlag', false);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 dim = 2;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 C_default = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 for l = 1:dim
1272
15865fbda16e Set default stiffness tensors to 0 in CurvilinearAnisotropic schemes
Martin Almquist <malmquist@stanford.edu>
parents: 1214
diff changeset
82 C_default{i,j,k,l} = @(x,y) 0*x ;
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 default_arg('C', C_default);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 assert(isa(g, 'grid.Curvilinear'));
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 if isa(rho, 'function_handle')
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 rho = grid.evalOn(g, rho);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95 C_mat = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 if isa(C{i,j,k,l}, 'function_handle')
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 obj.C = C_mat;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 m = g.size();
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 m_tot = g.N();
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113 % 1D operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 m_u = m(1);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 m_v = m(2);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 ops_u = opSet{1}(m_u, {0, 1}, order);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 ops_v = opSet{2}(m_v, {0, 1}, order);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119 h_u = ops_u.h;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 h_v = ops_v.h;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 I_u = speye(m_u);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 I_v = speye(m_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 D1_u = ops_u.D1;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 H_u = ops_u.H;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 Hi_u = ops_u.HI;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128 e_l_u = ops_u.e_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129 e_r_u = ops_u.e_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 d1_l_u = ops_u.d1_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 d1_r_u = ops_u.d1_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 D1_v = ops_v.D1;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134 H_v = ops_v.H;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135 Hi_v = ops_v.HI;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 e_l_v = ops_v.e_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 e_r_v = ops_v.e_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 d1_l_v = ops_v.d1_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139 d1_r_v = ops_v.d1_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142 % Logical operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143 Du = kr(D1_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144 Dv = kr(I_u,D1_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 e_w = kr(e_l_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 e_e = kr(e_r_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148 e_s = kr(I_u,e_l_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149 e_n = kr(I_u,e_r_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
150
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
151 % Metric coefficients
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
152 coords = g.points();
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 x = coords(:,1);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154 y = coords(:,2);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
156 x_u = Du*x;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
157 x_v = Dv*x;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
158 y_u = Du*y;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 y_v = Dv*y;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
160
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
161 J = x_u.*y_v - x_v.*y_u;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
162
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
163 K = cell(dim, dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
164 K{1,1} = y_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
165 K{1,2} = -y_u./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
166 K{2,1} = -x_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
167 K{2,2} = x_u./J;
1306
633757e582e5 Add transformation gradient as property in CurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1302
diff changeset
168 obj.K = K;
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
169
1211
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
170 % Physical derivatives
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
171 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
172 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
173
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174 % Wrap around Aniosotropic Cartesian
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 rho_tilde = J.*rho;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
176
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
177 PHI = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
182 PHI{i,j,k,l} = 0*C{i,j,k,l};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
183 for m = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184 for n = 1:dim
1209
67eee83fd9c9 Swap indices in anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1205
diff changeset
185 PHI{i,j,k,l} = PHI{i,j,k,l} + J.*K{m,i}.*C{m,j,n,l}.*K{n,k};
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
188 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
189 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
193 gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1});
1302
a0d615bde7f8 Add the hollow option to the anisotropic diffops
Martin Almquist <malmquist@stanford.edu>
parents: 1295
diff changeset
194 refObj = scheme.Elastic2dVariableAnisotropic(gRef, order, rho_tilde, PHI, opSet, [], hollow);
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196 %---- Set object properties ------
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
197 obj.RHO = spdiag(rho);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 % Volume quadrature
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200 obj.J = spdiag(J);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201 obj.Ji = spdiag(1./J);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 obj.H = obj.J*kr(H_u,H_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 % Boundary quadratures
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 obj.s_w = s_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210 obj.s_e = s_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 obj.s_s = s_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212 obj.s_n = s_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
213
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214 obj.H_w = H_v*spdiag(s_w);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 obj.H_e = H_v*spdiag(s_e);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 obj.H_s = H_u*spdiag(s_s);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 obj.H_n = H_u*spdiag(s_n);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
218
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219 % Restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 obj.e_w = refObj.e_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221 obj.e_e = refObj.e_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 obj.e_s = refObj.e_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 obj.e_n = refObj.e_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 % Adapt things from reference object
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226 obj.D = refObj.D;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
227 obj.E = refObj.E;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 obj.e1_w = refObj.e1_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 obj.e1_e = refObj.e1_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231 obj.e1_s = refObj.e1_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
232 obj.e1_n = refObj.e1_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
233
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
234 obj.e2_w = refObj.e2_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
235 obj.e2_e = refObj.e2_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
236 obj.e2_s = refObj.e2_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
237 obj.e2_n = refObj.e2_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
238
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
239 obj.e_scalar_w = refObj.e_scalar_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 obj.e_scalar_e = refObj.e_scalar_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241 obj.e_scalar_s = refObj.e_scalar_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
242 obj.e_scalar_n = refObj.e_scalar_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
243
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
244 e1_w = obj.e1_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
245 e1_e = obj.e1_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
246 e1_s = obj.e1_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
247 e1_n = obj.e1_n;
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
249 e2_w = obj.e2_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
250 e2_e = obj.e2_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
251 e2_s = obj.e2_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
252 e2_n = obj.e2_n;
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
254 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
257 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
258
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
259 obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
260 obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
261 obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
262 obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
263
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
264 obj.tau_w = (refObj.e_w'*obj.e1_w*obj.tau1_w')' + (refObj.e_w'*obj.e2_w*obj.tau2_w')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
265 obj.tau_e = (refObj.e_e'*obj.e1_e*obj.tau1_e')' + (refObj.e_e'*obj.e2_e*obj.tau2_e')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
266 obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
267 obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
268
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
269 % Physical normals
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
270 e_w = obj.e_scalar_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
271 e_e = obj.e_scalar_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
272 e_s = obj.e_scalar_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
273 e_n = obj.e_scalar_n;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
274
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
275 e_w_vec = obj.e_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
276 e_e_vec = obj.e_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
277 e_s_vec = obj.e_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
278 e_n_vec = obj.e_n;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
279
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
280 nu_w = [-1,0];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
281 nu_e = [1,0];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
282 nu_s = [0,-1];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
283 nu_n = [0,1];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
284
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
285 obj.n_w = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
286 obj.n_e = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
287 obj.n_s = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
288 obj.n_n = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
289
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
290 % Compute normal and rotate (exactly!) 90 degrees counter-clockwise to get tangent
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
291 n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
292 n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
293 obj.n_w{1} = spdiag(n_w_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
294 obj.n_w{2} = spdiag(n_w_2);
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
295 obj.tangent_w = {-obj.n_w{2}, obj.n_w{1}};
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
296
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
297 n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
298 n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
299 obj.n_e{1} = spdiag(n_e_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
300 obj.n_e{2} = spdiag(n_e_2);
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
301 obj.tangent_e = {-obj.n_e{2}, obj.n_e{1}};
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
302
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
303 n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
304 n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
305 obj.n_s{1} = spdiag(n_s_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
306 obj.n_s{2} = spdiag(n_s_2);
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
307 obj.tangent_s = {-obj.n_s{2}, obj.n_s{1}};
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
308
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
309 n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
310 n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
311 obj.n_n{1} = spdiag(n_n_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
312 obj.n_n{2} = spdiag(n_n_2);
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
313 obj.tangent_n = {-obj.n_n{2}, obj.n_n{1}};
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
314
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
315 % Operators that extract the normal component
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
316 obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
317 obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')';
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
318 obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
319 obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
320
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
321 % Operators that extract the tangential component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
322 obj.et_w = (obj.tangent_w{1}*obj.e1_w' + obj.tangent_w{2}*obj.e2_w')';
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
323 obj.et_e = (obj.tangent_e{1}*obj.e1_e' + obj.tangent_e{2}*obj.e2_e')';
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
324 obj.et_s = (obj.tangent_s{1}*obj.e1_s' + obj.tangent_s{2}*obj.e2_s')';
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
325 obj.et_n = (obj.tangent_n{1}*obj.e1_n' + obj.tangent_n{2}*obj.e2_n')';
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
326
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
327 obj.tau_n_w = (obj.n_w{1}*obj.tau1_w' + obj.n_w{2}*obj.tau2_w')';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
328 obj.tau_n_e = (obj.n_e{1}*obj.tau1_e' + obj.n_e{2}*obj.tau2_e')';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
329 obj.tau_n_s = (obj.n_s{1}*obj.tau1_s' + obj.n_s{2}*obj.tau2_s')';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
330 obj.tau_n_n = (obj.n_n{1}*obj.tau1_n' + obj.n_n{2}*obj.tau2_n')';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
331
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
332 obj.tau_t_w = (obj.tangent_w{1}*obj.tau1_w' + obj.tangent_w{2}*obj.tau2_w')';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
333 obj.tau_t_e = (obj.tangent_e{1}*obj.tau1_e' + obj.tangent_e{2}*obj.tau2_e')';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
334 obj.tau_t_s = (obj.tangent_s{1}*obj.tau1_s' + obj.tangent_s{2}*obj.tau2_s')';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
335 obj.tau_t_n = (obj.tangent_n{1}*obj.tau1_n' + obj.tangent_n{2}*obj.tau2_n')';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
336
1214
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
337 % Stress operators
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
338 sigma = cell(dim, dim);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
339 D1 = {obj.Dx, obj.Dy};
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
340 E = obj.E;
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
341 N = length(obj.RHO);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
342 for i = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
343 for j = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
344 sigma{i,j} = sparse(N,2*N);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
345 for k = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
346 for l = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
347 sigma{i,j} = sigma{i,j} + obj.C{i,j,k,l}*D1{k}*E{l}';
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
348 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
349 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
350 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
351 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
352 obj.sigma = sigma;
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
353
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
354 % Misc.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
355 obj.refObj = refObj;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
356 obj.m = refObj.m;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
357 obj.h = refObj.h;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
358 obj.order = order;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
359 obj.grid = g;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
360 obj.dim = dim;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
361
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
362 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
363
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
364
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
365 % Closure functions return the operators applied to the own domain to close the boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
366 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
367 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
368 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
369 % on the first component. Can also be e.g.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
370 % {'normal', 'd'} or {'tangential', 't'} for conditions on
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
371 % tangential/normal component.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
372 % data is a function returning the data that should be applied at the boundary.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
373 % neighbour_scheme is an instance of Scheme that should be interfaced to.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
374 % neighbour_boundary is a string specifying which boundary to interface to.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
375
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
376 % For displacement bc:
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377 % bc = {comp, 'd', dComps},
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
378 % where
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
379 % dComps = vector of components with displacement BC. Default: 1:dim.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
380 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
381 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
382 default_arg('tuning', 1.0);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
383 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
384
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
385 component = bc{1};
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
386 type = bc{2};
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
387
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
388 switch component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
389
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
390 % If conditions on Cartesian components
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
391 case {1,2}
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
392 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
393
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
394 % If conditions on normal or tangential components
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
395 case {'n', 't'}
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
396
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
397 switch component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
398 case 'n'
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
399 en = obj.getBoundaryOperator('en', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
400 case 't'
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
401 en = obj.getBoundaryOperator('et', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
402 end
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
403 e1 = obj.getBoundaryOperator('e1', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
404
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
405 bc1 = {1, type};
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
406 [c1, p1] = obj.refObj.boundary_condition(boundary, bc1, tuning);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
407 bc2 = {2, type};
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
408 c2 = obj.refObj.boundary_condition(boundary, bc2, tuning);
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
409
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
410 switch type
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
411 case {'F','f','Free','free','traction','Traction','t','T'}
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
412 closure = en*en'*(c1+c2);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
413 penalty = en*e1'*p1;
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
414 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
415 [closure, penalty] = obj.displacementBCNormalTangential(boundary, bc, tuning);
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
416 end
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
417
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
418 end
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
419
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
420 switch type
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
421 case {'F','f','Free','free','traction','Traction','t','T'}
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
422
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
423 s = obj.(['s_' boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
424 s = spdiag(s);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
425 penalty = penalty*s;
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
426
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
427 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
428 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
429
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
430 function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning)
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
431 disp('WARNING: displacementBCNormalTangential is only guaranteed to work for displaacement BC on one component and traction on the other');
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
432 u = obj;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
433
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
434 component = bc{1};
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
435 type = bc{2};
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
436
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
437 switch component
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
438 case 'n'
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
439 en = u.getBoundaryOperator('en', boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
440 tau_n = u.getBoundaryOperator('tau_n', boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
441 N = u.getNormal(boundary);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
442 case 't'
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
443 en = u.getBoundaryOperator('et', boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
444 tau_n = u.getBoundaryOperator('tau_t', boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
445 N = u.getTangent(boundary);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
446 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
447
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
448 % Operators
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
449 e = u.getBoundaryOperatorForScalarField('e', boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
450 h11 = u.getBorrowing(boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
451 n = u.getNormal(boundary);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
452
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
453 C = u.C;
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
454 Ji = u.Ji;
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
455 s = spdiag(u.(['s_' boundary]));
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
456 m_tot = u.grid.N();
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
457
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
458 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
459 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
460
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
461 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
462 dim = u.dim;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
463
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
464 % Preallocate
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
465 [~, m_int] = size(H_gamma);
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
466 closure = sparse(dim*m_tot, dim*m_tot);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
467 penalty = sparse(dim*m_tot, m_int);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
468
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
469 % Term 1: The symmetric term
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
470 Z = sparse(m_int, m_int);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
471 for i = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
472 for j = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
473 for l = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
474 for k = 1:dim
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
475 Z = Z + n{i}*N{j}*e'*Ji*C{i,j,k,l}*e*n{k}*N{l};
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
476 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
477 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
478 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
479 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
480
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
481 Z = -tuning*dim*1/h11*s*Z;
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
482 closure = closure + en*H_gamma*Z*en';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
483 penalty = penalty - en*H_gamma*Z;
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
484
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
485 % Term 2: The symmetrizing term
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
486 closure = closure + tau_n*H_gamma*en';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
487 penalty = penalty - tau_n*H_gamma;
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
488
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
489 % Multiply all normal component terms by inverse of density x quadraure
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
490 closure = RHOi*Hi*closure;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
491 penalty = RHOi*Hi*penalty;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
492 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
493
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
494 % type Struct that specifies the interface coupling.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
495 % Fields:
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
496 % -- tuning: penalty strength, defaults to 1.0
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
497 % -- interpolation: type of interpolation, default 'none'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
498 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
499
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
500 defaultType.tuning = 1.0;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
501 defaultType.interpolation = 'none';
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
502 defaultType.type = 'standard';
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
503 default_struct('type', defaultType);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
504
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
505 switch type.type
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
506 case 'standard'
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
507 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
508 case 'frictionalFault'
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
509 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
510 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
511
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
512 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
513
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
514 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
515 tuning = type.tuning;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
516
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
517 % u denotes the solution in the own domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
518 % v denotes the solution in the neighbour domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
519
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
520 u = obj;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
521 v = neighbour_scheme;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
522
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
523 % Operators, u side
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
524 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
525 en_u = u.getBoundaryOperator('en', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
526 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
527 h11_u = u.getBorrowing(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
528 n_u = u.getNormal(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
529
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
530 C_u = u.C;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
531 Ji_u = u.Ji;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
532 s_u = spdiag(u.(['s_' boundary]));
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
533 m_tot_u = u.grid.N();
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
534
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
535 % Operators, v side
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
536 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
537 en_v = v.getBoundaryOperator('en', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
538 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
539 h11_v = v.getBorrowing(neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
540 n_v = v.getNormal(neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
541
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
542 C_v = v.C;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
543 Ji_v = v.Ji;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
544 s_v = spdiag(v.(['s_' neighbour_boundary]));
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
545 m_tot_v = v.grid.N();
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
546
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
547 % Operators that are only required for own domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
548 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
549 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
550
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
551 % Shared operators
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
552 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
553 dim = u.dim;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
554
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
555 % Preallocate
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
556 [~, m_int] = size(H_gamma);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
557 closure = sparse(dim*m_tot_u, dim*m_tot_u);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
558 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
559
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
560 % Continuity of normal displacement, term 1: The symmetric term
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
561 Z_u = sparse(m_int, m_int);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
562 Z_v = sparse(m_int, m_int);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
563 for i = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
564 for j = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
565 for l = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
566 for k = 1:dim
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
567 Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*n_u{l};
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
568 Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*n_v{l};
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
569 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
570 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
571 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
572 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
573
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
574 Z = -tuning*dim*( 1/(4*h11_u)*s_u*Z_u + 1/(4*h11_v)*s_v*Z_v );
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
575 closure = closure + en_u*H_gamma*Z*en_u';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
576 penalty = penalty + en_u*H_gamma*Z*en_v';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
577
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
578 % Continuity of normal displacement, term 2: The symmetrizing term
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
579 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
580 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
581
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
582 % Continuity of normal traction
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
583 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
584 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
585
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
586 % Multiply all normal component terms by inverse of density x quadraure
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
587 closure = RHOi*Hi*closure;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
588 penalty = RHOi*Hi*penalty;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
589
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
590 % ---- Tangential tractions are imposed just like traction BC ------
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
591 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
592
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
593 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
594
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
595
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
596 % Returns h11 for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
597 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
598 function h11 = getBorrowing(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
599 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
600
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
601 switch boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
602 case {'w','e'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
603 h11 = obj.refObj.h11{1};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
604 case {'s', 'n'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
605 h11 = obj.refObj.h11{2};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
606 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
607 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
608
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
609 % Returns the outward unit normal vector for the boundary specified by the string boundary.
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
610 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
611 function n = getNormal(obj, boundary)
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
612 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
613
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
614 n = obj.(['n_' boundary]);
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
615 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
616
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
617 % Returns the unit tangent vector for the boundary specified by the string boundary.
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
618 % t is a cell of diagonal matrices for each normal component, t{1} = t_1, t{2} = t_2.
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
619 function t = getTangent(obj, boundary)
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
620 assertIsMember(boundary, {'w', 'e', 's', 'n'})
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
621
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
622 t = obj.(['tangent_' boundary]);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
623 end
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
624
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
625 % Returns the boundary operator op for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
626 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
627 function o = getBoundaryOperator(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
628 assertIsMember(boundary, {'w', 'e', 's', 'n'})
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
629 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
630
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
631 o = obj.([op, '_', boundary]);
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
632
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
633 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
634
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
635 % Returns the boundary operator op for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
636 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
637 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
638 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
639 assertIsMember(op, {'e'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
640
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
641 switch op
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
642
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
643 case 'e'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
644 o = obj.(['e_scalar', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
645 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
646
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
647 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
648
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
649 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
650 % Formula: tau_i = T_ij u_j
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
651 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
652 function T = getBoundaryTractionOperator(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
653 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
654
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
655 T = obj.(['T', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
656 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
657
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
658 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
659 % corresponding to the number of boundary unknowns
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
660 %
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
661 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
662 function H = getBoundaryQuadrature(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
663 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
664
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
665 H = obj.getBoundaryQuadratureForScalarField(boundary);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
666 I_dim = speye(obj.dim, obj.dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
667 H = kron(H, I_dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
668 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
669
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
670 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
671 % corresponding to the number of boundary grid points
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
672 %
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
673 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
674 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
675 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
676
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
677 H_b = obj.(['H_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
678 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
679
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
680 function N = size(obj)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
681 N = obj.dim*prod(obj.m);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
682 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
683 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
684 end