Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1302:a0d615bde7f8 feature/poroelastic

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