Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1214:e1d4cb8b5309 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add stress operators to scheme anisotropic curvilinear
author Martin Almquist <malmquist@stanford.edu>
date Mon, 21 Oct 2019 18:08:22 -0700
parents 43f1cd11e8e8
children 15865fbda16e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 classdef Elastic2dCurvilinearAnisotropic < scheme.Scheme
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 % Discretizes the elastic wave equation:
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5 % in curvilinear coordinates.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 % opSet should be cell array of opSets, one per dimension. This
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 % is useful if we have periodic BC in one direction.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 % Assumes fully compatible operators.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 properties
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 m % Number of points in each direction, possibly a vector
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12 h % Grid spacing
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14 grid
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 order % Order of accuracy for the approximation
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 % Diagonal matrices for variable coefficients
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 J, Ji
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21 RHO % Density
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 C % Elastic stiffness tensor
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24 D % Total operator
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25
1211
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
26 Dx, Dy % Physical derivatives
1214
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
27 sigma % Cell matrix of physical stress operators
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
28 n_w, n_e, n_s, n_n % Physical normals
1211
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
29
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 % Boundary operators in cell format, used for BC
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31 T_w, T_e, T_s, T_n
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 % Traction operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 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
35 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
36 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 % Inner products
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39 H
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 % Boundary inner products (for scalar field)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 H_w, H_e, H_s, H_n
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 % Surface Jacobian vectors
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45 s_w, s_e, s_s, s_n
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 % Boundary restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 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
49 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
50 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
51 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
52 en_w, en_e, en_s, en_n % Act on vector field, return normal component
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 % E{i}^T picks out component i
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55 E
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 % Elastic2dVariableAnisotropic object for reference domain
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 refObj
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 methods
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 % The coefficients can either be function handles or grid functions
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65 function obj = Elastic2dCurvilinearAnisotropic(g, order, rho, C, opSet, optFlag)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 default_arg('rho', @(x,y) 0*x+1);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 default_arg('optFlag', false);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 dim = 2;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 C_default = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 C_default{i,j,k,l} = @(x,y) 0*x + 1;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 default_arg('C', C_default);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 assert(isa(g, 'grid.Curvilinear'));
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 if isa(rho, 'function_handle')
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 rho = grid.evalOn(g, rho);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 end
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 C_mat = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 if isa(C{i,j,k,l}, 'function_handle')
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95 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
96 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 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
98 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 obj.C = C_mat;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 m = g.size();
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 m_tot = g.N();
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 % 1D operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 m_u = m(1);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 m_v = m(2);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 ops_u = opSet{1}(m_u, {0, 1}, order);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 ops_v = opSet{2}(m_v, {0, 1}, order);
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 h_u = ops_u.h;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 h_v = ops_v.h;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 I_u = speye(m_u);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 I_v = speye(m_v);
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 D1_u = ops_u.D1;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 H_u = ops_u.H;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 Hi_u = ops_u.HI;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 e_l_u = ops_u.e_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 e_r_u = ops_u.e_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 d1_l_u = ops_u.d1_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 d1_r_u = ops_u.d1_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 D1_v = ops_v.D1;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128 H_v = ops_v.H;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129 Hi_v = ops_v.HI;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 e_l_v = ops_v.e_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 e_r_v = ops_v.e_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132 d1_l_v = ops_v.d1_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 d1_r_v = ops_v.d1_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 % Logical operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 Du = kr(D1_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 Dv = kr(I_u,D1_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140 e_w = kr(e_l_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141 e_e = kr(e_r_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142 e_s = kr(I_u,e_l_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143 e_n = kr(I_u,e_r_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145 % Metric coefficients
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 coords = g.points();
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 x = coords(:,1);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148 y = coords(:,2);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
150 x_u = Du*x;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
151 x_v = Dv*x;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
152 y_u = Du*y;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 y_v = Dv*y;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155 J = x_u.*y_v - x_v.*y_u;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
156
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
157 K = cell(dim, dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
158 K{1,1} = y_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 K{1,2} = -y_u./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
160 K{2,1} = -x_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
161 K{2,2} = x_u./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
162
1211
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
163 % Physical derivatives
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic
Martin Almquist <malmquist@stanford.edu>
parents: 1209
diff changeset
164 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
165 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
166
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
167 % Wrap around Aniosotropic Cartesian
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168 rho_tilde = J.*rho;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
169
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
170 PHI = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
171 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
172 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 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
176 for m = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
177 for n = 1:dim
1209
67eee83fd9c9 Swap indices in anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1205
diff changeset
178 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
179 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
182 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
183 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
185
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186 gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1});
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187 refObj = scheme.Elastic2dVariableAnisotropic(gRef, order, rho_tilde, PHI, opSet);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
188
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
189 %---- Set object properties ------
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 obj.RHO = spdiag(rho);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192 % Volume quadrature
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
193 obj.J = spdiag(J);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
194 obj.Ji = spdiag(1./J);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195 obj.H = obj.J*kr(H_u,H_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
197 % Boundary quadratures
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198 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
199 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
200 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
201 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
202 obj.s_w = s_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 obj.s_e = s_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 obj.s_s = s_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 obj.s_n = s_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207 obj.H_w = H_v*spdiag(s_w);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 obj.H_e = H_v*spdiag(s_e);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 obj.H_s = H_u*spdiag(s_s);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210 obj.H_n = H_u*spdiag(s_n);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212 % Restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
213 obj.e_w = refObj.e_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214 obj.e_e = refObj.e_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 obj.e_s = refObj.e_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 obj.e_n = refObj.e_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 % Adapt things from reference object
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219 obj.D = refObj.D;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 obj.E = refObj.E;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 obj.e1_w = refObj.e1_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 obj.e1_e = refObj.e1_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 obj.e1_s = refObj.e1_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 obj.e1_n = refObj.e1_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
227 obj.e2_w = refObj.e2_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228 obj.e2_e = refObj.e2_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 obj.e2_s = refObj.e2_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 obj.e2_n = refObj.e2_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
232 obj.e_scalar_w = refObj.e_scalar_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
233 obj.e_scalar_e = refObj.e_scalar_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
234 obj.e_scalar_s = refObj.e_scalar_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
235 obj.e_scalar_n = refObj.e_scalar_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
236
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
237 e1_w = obj.e1_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
238 e1_e = obj.e1_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
239 e1_s = obj.e1_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
240 e1_n = obj.e1_n;
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
242 e2_w = obj.e2_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
243 e2_e = obj.e2_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
244 e2_s = obj.e2_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
245 e2_n = obj.e2_n;
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
246
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
247 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
249 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
250 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
251
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
252 obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253 obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
254 obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
257 obj.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
258 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
259 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
260 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
261
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
262 % Physical normals
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
263 e_w = obj.e_scalar_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
264 e_e = obj.e_scalar_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
265 e_s = obj.e_scalar_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
266 e_n = obj.e_scalar_n;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
267
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
268 e_w_vec = obj.e_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
269 e_e_vec = obj.e_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
270 e_s_vec = obj.e_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
271 e_n_vec = obj.e_n;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
272
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
273 nu_w = [-1,0];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
274 nu_e = [1,0];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
275 nu_s = [0,-1];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
276 nu_n = [0,1];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
277
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
278 obj.n_w = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
279 obj.n_e = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
280 obj.n_s = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
281 obj.n_n = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
282
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
283 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
284 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
285 obj.n_w{1} = spdiag(n_w_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
286 obj.n_w{2} = spdiag(n_w_2);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
287
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
288 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
289 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
290 obj.n_e{1} = spdiag(n_e_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
291 obj.n_e{2} = spdiag(n_e_2);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
292
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
293 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
294 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
295 obj.n_s{1} = spdiag(n_s_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
296 obj.n_s{2} = spdiag(n_s_2);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
297
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
298 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
299 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
300 obj.n_n{1} = spdiag(n_n_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
301 obj.n_n{2} = spdiag(n_n_2);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
302
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
303 % Operators that extract the normal component
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
304 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
305 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
306 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
307 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
308
1214
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
309 % Stress operators
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
310 sigma = cell(dim, dim);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
311 D1 = {obj.Dx, obj.Dy};
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
312 E = obj.E;
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
313 N = length(obj.RHO);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
314 for i = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
315 for j = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
316 sigma{i,j} = sparse(N,2*N);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
317 for k = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
318 for l = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
319 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
320 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
321 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
322 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
323 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
324 obj.sigma = sigma;
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1213
diff changeset
325
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
326 % Misc.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
327 obj.refObj = refObj;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
328 obj.m = refObj.m;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
329 obj.h = refObj.h;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
330 obj.order = order;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
331 obj.grid = g;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
332 obj.dim = dim;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
333
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
334 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
335
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
336
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
337 % 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
338 % 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
339 % 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
340 % 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
341 % on the first component. Can also be e.g.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
342 % {'normal', 'd'} or {'tangential', 't'} for conditions on
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
343 % tangential/normal component.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
344 % 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
345 % 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
346 % 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
347
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
348 % For displacement bc:
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
349 % bc = {comp, 'd', dComps},
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
350 % where
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 % dComps = vector of components with displacement BC. Default: 1:dim.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352 % 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
353 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
354 default_arg('tuning', 1.0);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
355 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
356
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
357 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
358
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
359 type = bc{2};
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 switch type
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
362 case {'F','f','Free','free','traction','Traction','t','T'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
363 s = obj.(['s_' boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
364 s = spdiag(s);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
365 penalty = penalty*s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
366 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
367 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
368
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
369 % type Struct that specifies the interface coupling.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
370 % Fields:
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
371 % -- tuning: penalty strength, defaults to 1.0
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
372 % -- interpolation: type of interpolation, default 'none'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
373 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
374
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
375 defaultType.tuning = 1.0;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
376 defaultType.interpolation = 'none';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377 default_struct('type', defaultType);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
378
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
379 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
380 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
381
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
382 % Returns h11 for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
383 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
384 function h11 = getBorrowing(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
385 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
386
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
387 switch boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
388 case {'w','e'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
389 h11 = obj.refObj.h11{1};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
390 case {'s', 'n'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
391 h11 = obj.refObj.h11{2};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
392 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
393 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
394
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
395 % 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
396 % 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
397 function n = getNormal(obj, boundary)
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
398 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
399
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
400 n = obj.(['n_' boundary]);
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
401 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
402
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
403 % 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
404 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
405 function o = getBoundaryOperator(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
406 assertIsMember(boundary, {'w', 'e', 's', 'n'})
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
407 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'})
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
408
1213
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents: 1211
diff changeset
409 o = obj.([op, '_', boundary]);
1205
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
410
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
411 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
412
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
413 % 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
414 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
415 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
416 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
417 assertIsMember(op, {'e'})
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 op
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
420
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
421 case 'e'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
422 o = obj.(['e_scalar', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
423 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
424
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
425 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
426
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
427 % 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
428 % Formula: tau_i = T_ij u_j
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
429 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
430 function T = getBoundaryTractionOperator(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
431 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
432
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
433 T = obj.(['T', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
434 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
435
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
436 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
437 % corresponding to the number of boundary unknowns
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
438 %
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
439 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
440 function H = getBoundaryQuadrature(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
441 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
442
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
443 H = obj.getBoundaryQuadratureForScalarField(boundary);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
444 I_dim = speye(obj.dim, obj.dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
445 H = kron(H, I_dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
446 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
447
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
448 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
449 % corresponding to the number of boundary grid points
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
450 %
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
451 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
452 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
453 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
454
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
455 H_b = obj.(['H_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
456 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
457
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
458 function N = size(obj)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
459 N = obj.dim*prod(obj.m);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
460 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
461 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
462 end