Mercurial > repos > public > sbplib

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

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