Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1338:da61892884a4 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for using periodic grids.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 13 May 2022 13:29:05 +0200
parents df8c71b80c33
children
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);
1338
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
155
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
156 % When computing the metric derivatives, we can't use periodic operators
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
157 if isequal(opSet{1},@sbp.D2VariablePeriodic)
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
158 ops_metric = sbp.D2Standard(m_u, {0, 1-1/m_u}, order);
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
159 Du_metric = kr(ops_metric.D1,I_v);
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
160 else
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
161 Du_metric = Du;
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
162 end
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
163
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
164 if isequal(opSet{2},@sbp.D2VariablePeriodic)
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
165 ops_metric = sbp.D2Standard(m_v, {0, 1-1/m_v}, order);
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
166 Dv_metric = kr(I_u,ops_metric.D1);
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
167 else
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
168 Dv_metric = Dv;
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
169 end
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
170
1338
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
171 x_u = Du_metric*x;
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
172 x_v = Dv_metric*x;
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
173 y_u = Du_metric*y;
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
174 y_v = Dv_metric*y;
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
176 J = x_u.*y_v - x_v.*y_u;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
177
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178 K = cell(dim, dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179 K{1,1} = y_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 K{1,2} = -y_u./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181 K{2,1} = -x_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
182 K{2,2} = x_u./J;
1306
633757e582e5 Add transformation gradient as property in CurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1302
diff changeset
183 obj.K = K;
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184
1211
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
185 % Physical derivatives
1338
da61892884a4 Add support for using periodic grids.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1334
diff changeset
186 % TODO: Fix for periodic operators
1211
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
187 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
188 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
189
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 % Wrap around Aniosotropic Cartesian
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191 rho_tilde = J.*rho;
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 PHI = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
194 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
197 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198 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
199 for m = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200 for n = 1:dim
1209
67eee83fd9c9 Swap indices in anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1205
diff changeset
201 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
202 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208
1334
df8c71b80c33 Use the logic grid associated with a CurvilinearGrid instead of creating a new one
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1322
diff changeset
209 refObj = scheme.Elastic2dVariableAnisotropic(g.logic, order, rho_tilde, PHI, opSet, [], hollow);
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 %---- Set object properties ------
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212 obj.RHO = spdiag(rho);
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 % Volume quadrature
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 obj.J = spdiag(J);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 obj.Ji = spdiag(1./J);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 obj.H = obj.J*kr(H_u,H_v);
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 % Boundary quadratures
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 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
221 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
222 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
223 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
224 obj.s_w = s_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 obj.s_e = s_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226 obj.s_s = s_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
227 obj.s_n = s_n;
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.H_w = H_v*spdiag(s_w);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 obj.H_e = H_v*spdiag(s_e);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231 obj.H_s = H_u*spdiag(s_s);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
232 obj.H_n = H_u*spdiag(s_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 % Restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
235 obj.e_w = refObj.e_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
236 obj.e_e = refObj.e_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
237 obj.e_s = refObj.e_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
238 obj.e_n = refObj.e_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
239
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 % Adapt things from reference object
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241 obj.D = refObj.D;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
242 obj.E = refObj.E;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
243
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
244 obj.e1_w = refObj.e1_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
245 obj.e1_e = refObj.e1_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
246 obj.e1_s = refObj.e1_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
247 obj.e1_n = refObj.e1_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
249 obj.e2_w = refObj.e2_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
250 obj.e2_e = refObj.e2_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
251 obj.e2_s = refObj.e2_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
252 obj.e2_n = refObj.e2_n;
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.e_scalar_w = refObj.e_scalar_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 obj.e_scalar_e = refObj.e_scalar_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256 obj.e_scalar_s = refObj.e_scalar_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
257 obj.e_scalar_n = refObj.e_scalar_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
258
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
259 e1_w = obj.e1_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
260 e1_e = obj.e1_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
261 e1_s = obj.e1_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
262 e1_n = obj.e1_n;
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
263
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
264 e2_w = obj.e2_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
265 e2_e = obj.e2_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
266 e2_s = obj.e2_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
267 e2_n = obj.e2_n;
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
268
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
269 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
271 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
272 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
273
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
274 obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
275 obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
276 obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
277 obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
278
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
279 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
280 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
281 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
282 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
283
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
284 % Physical normals
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
285 e_w = obj.e_scalar_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
286 e_e = obj.e_scalar_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
287 e_s = obj.e_scalar_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
288 e_n = obj.e_scalar_n;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
289
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
290 e_w_vec = obj.e_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
291 e_e_vec = obj.e_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
292 e_s_vec = obj.e_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
293 e_n_vec = obj.e_n;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
294
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
295 nu_w = [-1,0];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
296 nu_e = [1,0];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
297 nu_s = [0,-1];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
298 nu_n = [0,1];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
299
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
300 obj.n_w = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
301 obj.n_e = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
302 obj.n_s = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
303 obj.n_n = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
304
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
305 % 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
306 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
307 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
308 obj.n_w{1} = spdiag(n_w_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
309 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
310 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
311
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
312 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
313 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
314 obj.n_e{1} = spdiag(n_e_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
315 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
316 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
317
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
318 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
319 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
320 obj.n_s{1} = spdiag(n_s_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
321 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
322 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
323
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
324 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
325 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
326 obj.n_n{1} = spdiag(n_n_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
327 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
328 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
329
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
330 % Operators that extract the normal component
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
331 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
332 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
333 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
334 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
335
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
336 % Operators that extract the tangential component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
337 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
338 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
339 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
340 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
341
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
342 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
343 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
344 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
345 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
346
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
347 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
348 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
349 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
350 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
351
1214
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
352 % Stress operators
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
353 sigma = cell(dim, dim);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
354 D1 = {obj.Dx, obj.Dy};
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
355 E = obj.E;
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
356 N = length(obj.RHO);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
357 for i = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
358 for j = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
359 sigma{i,j} = sparse(N,2*N);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
360 for k = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
361 for l = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
362 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
363 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
364 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
365 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
366 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
367 obj.sigma = sigma;
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
368
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
369 % Misc.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
370 obj.refObj = refObj;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
371 obj.m = refObj.m;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
372 obj.h = refObj.h;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
373 obj.order = order;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
374 obj.grid = g;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
375 obj.dim = dim;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
376
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
378
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
379
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
380 % 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
381 % 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
382 % 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
383 % 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
384 % on the first component. Can also be e.g.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
385 % {'normal', 'd'} or {'tangential', 't'} for conditions on
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
386 % tangential/normal component.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
387 % 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
388 % 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
389 % 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
390
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
391 % For displacement bc:
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
392 % bc = {comp, 'd', dComps},
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
393 % where
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
394 % dComps = vector of components with displacement BC. Default: 1:dim.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
395 % 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
396 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
397 default_arg('tuning', 1.0);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
398 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
399
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
400 component = bc{1};
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
401 type = bc{2};
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
402
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
403 switch component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
404
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
405 % If conditions on Cartesian components
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
406 case {1,2}
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
407 [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
408
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
409 % If conditions on normal or tangential components
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
410 case {'n', 't'}
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
411
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
412 switch component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
413 case 'n'
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
414 en = obj.getBoundaryOperator('en', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
415 case 't'
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
416 en = obj.getBoundaryOperator('et', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
417 end
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
418 e1 = obj.getBoundaryOperator('e1', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
419
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
420 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
421 [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
422 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
423 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
424
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
425 switch type
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
426 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
427 closure = en*en'*(c1+c2);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
428 penalty = en*e1'*p1;
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
429 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
430 [closure, penalty] = obj.displacementBCNormalTangential(boundary, bc, tuning);
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
431 end
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
432
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
433 end
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
434
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
435 switch type
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
436 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
437
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
438 s = obj.(['s_' boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
439 s = spdiag(s);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
440 penalty = penalty*s;
1291
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
441
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
442 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
443 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
444
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
445 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
446 u = obj;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
447
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
448 component = bc{1};
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
449 type = bc{2};
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
450
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
451 switch component
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
452 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
453 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
454 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
455 N = u.getNormal(boundary);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
456 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
457 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
458 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
459 N = u.getTangent(boundary);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
460 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
461
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
462 % 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
463 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
464 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
465 n = u.getNormal(boundary);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
466
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
467 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
468 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
469 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
470 m_tot = u.grid.N();
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
471
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
472 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
473 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
474
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
475 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
476 dim = u.dim;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
477
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
478 % Preallocate
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
479 [~, 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
480 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
481 penalty = sparse(dim*m_tot, m_int);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
482
1295
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.
Martin Almquist <malmquist@stanford.edu>
parents: 1294
diff changeset
483 % Term 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
484 Z = sparse(m_int, m_int);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
485 for i = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
486 for j = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
487 for l = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
488 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
489 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
490 end
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 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
493 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
494
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
495 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
496 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
497 penalty = penalty - en*H_gamma*Z;
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
498
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
499 % 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
500 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
501 penalty = penalty - tau_n*H_gamma;
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
502
1315
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
503 % 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
504 closure = RHOi*Hi*closure;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
505 penalty = RHOi*Hi*penalty;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
506 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
507
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
508 % type Struct that specifies the interface coupling.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
509 % Fields:
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
510 % -- tuning: penalty strength, defaults to 1.0
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
511 % -- interpolation: type of interpolation, default 'none'
1317
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
512 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
513
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
514 defaultType.tuning = 1.0;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
515 defaultType.interpolation = 'none';
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
516 defaultType.type = 'standard';
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
517 default_struct('type', defaultType);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
518
1317
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
519 forcingPenalties = [];
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
520
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
521 switch type.type
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
522 case 'standard'
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
523 [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
524 case 'normalTangential'
1317
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
525 [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
526 case 'frictionalFault'
1317
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
527 [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
528 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
529
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
530 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
531
1317
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
532 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
533 tuning = type.tuning;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
534
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
535 % u denotes the solution in the own domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
536 % v denotes the solution in the neighbour domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
537
1317
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
538 forcingPenalties = cell(1, 1);
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
539 u = obj;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
540 v = neighbour_scheme;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
541
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
542 % Operators, u side
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
543 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
544 en_u = u.getBoundaryOperator('en', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
545 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
546 h11_u = u.getBorrowing(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
547 n_u = u.getNormal(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
548
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
549 C_u = u.C;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
550 Ji_u = u.Ji;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
551 s_u = spdiag(u.(['s_' boundary]));
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
552 m_tot_u = u.grid.N();
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
553
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
554 % Operators, v side
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
555 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
556 en_v = v.getBoundaryOperator('en', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
557 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
558 h11_v = v.getBorrowing(neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
559 n_v = v.getNormal(neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
560
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
561 C_v = v.C;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
562 Ji_v = v.Ji;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
563 s_v = spdiag(v.(['s_' neighbour_boundary]));
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
564 m_tot_v = v.grid.N();
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
565
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
566 % Operators that are only required for own domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
567 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
568 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
569
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
570 % Shared operators
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
571 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
572 dim = u.dim;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
573
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
574 % Preallocate
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
575 [~, m_int] = size(H_gamma);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
576 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
577 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
578
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
579 % 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
580 Z_u = sparse(m_int, m_int);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
581 Z_v = sparse(m_int, m_int);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
582 for i = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
583 for j = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
584 for l = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
585 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
586 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
587 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
588 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
589 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
590 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
591 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
592
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
593 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
594 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
595 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
596
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
597 % 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
598 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
599 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
600
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
601 % Continuity of normal traction
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
602 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
603 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
604 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
605
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
606 % 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
607 closure = RHOi*Hi*closure;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
608 penalty = RHOi*Hi*penalty;
1317
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
609 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
610
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
611 % ---- Tangential tractions are imposed just like traction BC ------
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
612 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
613
1317
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
614 forcingPenalties{1} = forcing_tau_n;
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
615
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
616 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
617
1315
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
618 % 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
619 % 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
620 % The resulting scheme is not identical to that of interfaceStandard. This appears to be better.
1322
412b8ceafbc6 Add Zt to output args for Elastic2dCurvilinearAnisotropic.interfaceNormalTangential
Martin Almquist <malmquist@stanford.edu>
parents: 1321
diff changeset
621 function [closure, penalty, forcingPenalties, Zt] = 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
622 tuning = type.tuning;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
623
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
624 % 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
625 % 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
626
1317
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
627 forcingPenalties = cell(2, 1);
1315
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
628 u = obj;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
629 v = neighbour_scheme;
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, u side
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
632 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
633 en_u = u.getBoundaryOperator('en', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
634 et_u = u.getBoundaryOperator('et', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
635 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
636 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
637 h11_u = u.getBorrowing(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
638 n_u = u.getNormal(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
639 t_u = u.getTangent(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_u = u.C;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
642 Ji_u = u.Ji;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
643 s_u = spdiag(u.(['s_' boundary]));
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
644 m_tot_u = u.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, v side
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
647 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
648 en_v = v.getBoundaryOperator('en', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
649 et_v = v.getBoundaryOperator('et', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
650 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
651 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
652 h11_v = v.getBorrowing(neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
653 n_v = v.getNormal(neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
654 t_v = v.getTangent(neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
655
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
656 C_v = v.C;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
657 Ji_v = v.Ji;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
658 s_v = spdiag(v.(['s_' neighbour_boundary]));
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
659 m_tot_v = v.grid.N();
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
660
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
661 % 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
662 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
663 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
664
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
665 % Shared operators
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
666 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
667 dim = u.dim;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
668
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
669 % Preallocate
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
670 [~, m_int] = size(H_gamma);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
671 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
672 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
673
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
674 % -- 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
675 Zn_u = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
676 Zn_v = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
677 Zt_u = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
678 Zt_v = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
679 for i = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
680 for j = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
681 for l = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
682 for k = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
683 % Penalty strength for normal component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
684 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
685 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
686
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
687 % Penalty strength for tangential component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
688 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
689 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
690 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
691 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
692 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
693 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
694
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
695 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
696 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
697
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
698 % Continuity of normal component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
699 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
700 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
701 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
702
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
703 % Continuity of tangential component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
704 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
705 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
706 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
707 %------------------------------------------------------------------
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 displacement, term 2: The symmetrizing term
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
710
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
711 % Continuity of normal displacement
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
712 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
713 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
714 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
715
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
716 % Continuity of tangential displacement
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
717 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
718 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
719 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
720 % ------------------------------------------------------------------
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
721
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
722 % --- Continuity of tractions -----------------------------
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
723
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
724 % Continuity of normal traction
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
725 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
726 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
727
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
728 % Continuity of tangential traction
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
729 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
730 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
731 %--------------------------------------------------------------------
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
732
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
733 % 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
734 closure = RHOi*Hi*closure;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
735 penalty = RHOi*Hi*penalty;
1317
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
736 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
737 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
738
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
739 forcingPenalties{1} = forcing_u_n;
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1315
diff changeset
740 forcingPenalties{2} = forcing_u_t;
1315
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
741
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
742 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
Martin Almquist <malmquist@stanford.edu>
parents: 1306
diff changeset
743
1294
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve
Martin Almquist <malmquist@stanford.edu>
parents: 1291
diff changeset
744
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
745 % Returns h11 for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
746 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
747 function h11 = getBorrowing(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
748 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
749
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
750 switch boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
751 case {'w','e'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
752 h11 = obj.refObj.h11{1};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
753 case {'s', 'n'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
754 h11 = obj.refObj.h11{2};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
755 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
756 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
757
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
758 % 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
759 % 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
760 function n = getNormal(obj, boundary)
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
761 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
762
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
763 n = obj.(['n_' boundary]);
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
764 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
765
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
766 % 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
767 % 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
768 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
769 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
770
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
771 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
772 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
773
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
774 % 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
775 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
776 function o = getBoundaryOperator(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
777 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
778 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
779
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
780 o = obj.([op, '_', boundary]);
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
781
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
782 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
783
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
784 % 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
785 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
786 function o = getBoundaryOperatorForScalarField(obj, op, 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 assertIsMember(op, {'e'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
789
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
790 switch op
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 case 'e'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
793 o = obj.(['e_scalar', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
794 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
795
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
796 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
797
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
798 % 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
799 % Formula: tau_i = T_ij u_j
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
800 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
801 function T = getBoundaryTractionOperator(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
802 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 T = obj.(['T', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
805 end
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 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
808 % corresponding to the number of boundary unknowns
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
809 %
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
810 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
811 function H = getBoundaryQuadrature(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
812 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 H = obj.getBoundaryQuadratureForScalarField(boundary);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
815 I_dim = speye(obj.dim, obj.dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
816 H = kron(H, I_dim);
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
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
819 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
820 % corresponding to the number of boundary grid points
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
821 %
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
822 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
823 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
824 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
825
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
826 H_b = obj.(['H_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
827 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
828
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
829 function N = size(obj)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
830 N = obj.dim*prod(obj.m);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
831 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
832 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
833 end