Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1321:b40359c9faed feature/poroelastic

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