Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1205:3258dca12af8 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add scheme for anisotropic curvilinear
author Martin Almquist <malmquist@stanford.edu>
date Sat, 07 Sep 2019 13:05:26 -0700
parents
children 67eee83fd9c9
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
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 % Boundary operators in cell format, used for BC
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 T_w, T_e, T_s, T_n
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29 % Traction operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 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
31 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
32 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
33
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 % Inner products
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35 H
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 % Boundary inner products (for scalar field)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 H_w, H_e, H_s, H_n
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 % Surface Jacobian vectors
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 s_w, s_e, s_s, s_n
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 % Boundary restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 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
45 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
46 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
47 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 % E{i}^T picks out component i
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 E
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 % Elastic2dVariableAnisotropic object for reference domain
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 refObj
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56 methods
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 % The coefficients can either be function handles or grid functions
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 % 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
60 function obj = Elastic2dCurvilinearAnisotropic(g, order, rho, C, opSet, optFlag)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 default_arg('rho', @(x,y) 0*x+1);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 default_arg('optFlag', false);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 dim = 2;
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 C_default = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 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
72 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 default_arg('C', C_default);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 assert(isa(g, 'grid.Curvilinear'));
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 if isa(rho, 'function_handle')
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 rho = grid.evalOn(g, rho);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 C_mat = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 if isa(C{i,j,k,l}, 'function_handle')
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 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
91 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 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
93 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95 end
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 obj.C = C_mat;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 m = g.size();
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 m_tot = g.N();
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 % 1D operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 m_u = m(1);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 m_v = m(2);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 ops_u = opSet{1}(m_u, {0, 1}, order);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 ops_v = opSet{2}(m_v, {0, 1}, order);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 h_u = ops_u.h;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 h_v = ops_v.h;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 I_u = speye(m_u);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 I_v = speye(m_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 D1_u = ops_u.D1;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 H_u = ops_u.H;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 Hi_u = ops_u.HI;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 e_l_u = ops_u.e_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118 e_r_u = ops_u.e_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119 d1_l_u = ops_u.d1_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 d1_r_u = ops_u.d1_r;
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 D1_v = ops_v.D1;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 H_v = ops_v.H;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 Hi_v = ops_v.HI;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 e_l_v = ops_v.e_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 e_r_v = ops_v.e_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 d1_l_v = ops_v.d1_l;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128 d1_r_v = ops_v.d1_r;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 % Logical operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132 Du = kr(D1_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 Dv = kr(I_u,D1_v);
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 e_w = kr(e_l_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 e_e = kr(e_r_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 e_s = kr(I_u,e_l_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 e_n = kr(I_u,e_r_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 % Metric coefficients
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141 coords = g.points();
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142 x = coords(:,1);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143 y = coords(:,2);
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 x_u = Du*x;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 x_v = Dv*x;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 y_u = Du*y;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148 y_v = Dv*y;
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 J = x_u.*y_v - x_v.*y_u;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
151
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
152 K = cell(dim, dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 K{1,1} = y_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154 K{1,2} = -y_u./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155 K{2,1} = -x_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
156 K{2,2} = x_u./J;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
157
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
158 % Wrap around Aniosotropic Cartesian
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 rho_tilde = J.*rho;
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 PHI = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
162 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
163 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
164 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
165 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
166 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
167 for m = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168 for n = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
169 PHI{i,j,k,l} = PHI{i,j,k,l} + J.*K{m,j}.*C{i,m,n,l}.*K{n,k};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
170 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
171 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
172 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 end
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 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
178 refObj = scheme.Elastic2dVariableAnisotropic(gRef, order, rho_tilde, PHI, opSet);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 %---- Set object properties ------
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181 obj.RHO = spdiag(rho);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
182
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
183 % Volume quadrature
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184 obj.J = spdiag(J);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
185 obj.Ji = spdiag(1./J);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186 obj.H = obj.J*kr(H_u,H_v);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
188 % Boundary quadratures
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
189 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
190 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
191 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
192 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
193 obj.s_w = s_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
194 obj.s_e = s_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195 obj.s_s = s_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196 obj.s_n = s_n;
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 obj.H_w = H_v*spdiag(s_w);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 obj.H_e = H_v*spdiag(s_e);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200 obj.H_s = H_u*spdiag(s_s);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201 obj.H_n = H_u*spdiag(s_n);
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 % Restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 obj.e_w = refObj.e_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 obj.e_e = refObj.e_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206 obj.e_s = refObj.e_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207 obj.e_n = refObj.e_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 % Adapt things from reference object
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210 obj.D = refObj.D;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 obj.E = refObj.E;
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.e1_w = refObj.e1_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214 obj.e1_e = refObj.e1_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 obj.e1_s = refObj.e1_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 obj.e1_n = refObj.e1_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 obj.e2_w = refObj.e2_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219 obj.e2_e = refObj.e2_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 obj.e2_s = refObj.e2_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221 obj.e2_n = refObj.e2_n;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 obj.e_scalar_w = refObj.e_scalar_w;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 obj.e_scalar_e = refObj.e_scalar_e;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 obj.e_scalar_s = refObj.e_scalar_s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226 obj.e_scalar_n = refObj.e_scalar_n;
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 e1_w = (obj.e_scalar_w'*obj.E{1}')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 e1_e = (obj.e_scalar_e'*obj.E{1}')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 e1_s = (obj.e_scalar_s'*obj.E{1}')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231 e1_n = (obj.e_scalar_n'*obj.E{1}')';
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 e2_w = (obj.e_scalar_w'*obj.E{2}')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
234 e2_e = (obj.e_scalar_e'*obj.E{2}')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
235 e2_s = (obj.e_scalar_s'*obj.E{2}')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
236 e2_n = (obj.e_scalar_n'*obj.E{2}')';
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.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
239 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
242
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
243 obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
244 obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
245 obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
246 obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
247
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248 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
249 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
250 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
251 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
252
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253 % Misc.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
254 obj.refObj = refObj;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 obj.m = refObj.m;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256 obj.h = refObj.h;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
257 obj.order = order;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
258 obj.grid = g;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
259 obj.dim = dim;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
260
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
261 end
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
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
264 % 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
265 % 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
266 % 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
267 % 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
268 % on the first component. Can also be e.g.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
269 % {'normal', 'd'} or {'tangential', 't'} for conditions on
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270 % tangential/normal component.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
271 % 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
272 % 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
273 % 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
274
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
275 % For displacement bc:
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
276 % bc = {comp, 'd', dComps},
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
277 % where
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
278 % dComps = vector of components with displacement BC. Default: 1:dim.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
279 % 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
280 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
281 default_arg('tuning', 1.0);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
282 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
283
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
284 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
285
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
286 type = bc{2};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
287
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
288 switch type
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
289 case {'F','f','Free','free','traction','Traction','t','T'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
290 s = obj.(['s_' boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
291 s = spdiag(s);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
292 penalty = penalty*s;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
293 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
294 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
295
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
296 % type Struct that specifies the interface coupling.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
297 % Fields:
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
298 % -- tuning: penalty strength, defaults to 1.0
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
299 % -- interpolation: type of interpolation, default 'none'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
300 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
301
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
302 defaultType.tuning = 1.0;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
303 defaultType.interpolation = 'none';
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
304 default_struct('type', defaultType);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
305
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
306 [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
307 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
308
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
309 % Returns the component number that is the tangential/normal component
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
310 % at the specified boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
311 function comp = getComponent(obj, comp_str, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
312 assertIsMember(comp_str, {'normal', 'tangential'});
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
313 assertIsMember(boundary, {'w', 'e', 's', 'n'});
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
314
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
315 switch boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
316 case {'w', 'e'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
317 switch comp_str
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
318 case 'normal'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
319 comp = 1;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
320 case 'tangential'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
321 comp = 2;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
322 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
323 case {'s', 'n'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
324 switch comp_str
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
325 case 'normal'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
326 comp = 2;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
327 case 'tangential'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
328 comp = 1;
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
329 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
330 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
331 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
332
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
333 % Returns h11 for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
334 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
335 function h11 = getBorrowing(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
336 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
337
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
338 switch boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
339 case {'w','e'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
340 h11 = obj.refObj.h11{1};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
341 case {'s', 'n'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
342 h11 = obj.refObj.h11{2};
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
343 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
344 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
345
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
346 % Returns the outward unit normal vector for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
347 function nu = getNormal(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
348 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
349
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
350 switch boundary
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 case 'w'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352 nu = [-1,0];
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
353 case 'e'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
354 nu = [1,0];
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
355 case 's'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
356 nu = [0,-1];
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
357 case 'n'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
358 nu = [0,1];
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
359 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
360 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
361
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
362 % 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
363 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
364 function o = getBoundaryOperator(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
365 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
366 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
367
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
368 switch op
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
369 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
370 o = obj.([op, '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
371 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
372
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
373 end
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 % 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
376 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
378 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
379 assertIsMember(op, {'e'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
380
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
381 switch op
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
382
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
383 case 'e'
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
384 o = obj.(['e_scalar', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
385 end
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 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
388
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
389 % 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
390 % Formula: tau_i = T_ij u_j
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
391 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
392 function T = getBoundaryTractionOperator(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
393 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 T = obj.(['T', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
396 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
397
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
398 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
399 % corresponding to the number of boundary unknowns
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
400 %
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
401 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
402 function H = getBoundaryQuadrature(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
403 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
404
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
405 H = obj.getBoundaryQuadratureForScalarField(boundary);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
406 I_dim = speye(obj.dim, obj.dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
407 H = kron(H, I_dim);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
408 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
409
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
410 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
411 % corresponding to the number of boundary grid points
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 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
414 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
415 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
416
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
417 H_b = obj.(['H_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
418 end
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
419
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
420 function N = size(obj)
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
421 N = obj.dim*prod(obj.m);
3258dca12af8 Add scheme for anisotropic curvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
422 end
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 end