Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dVariableAnisotropicMixedStencil.m @ 1339:bcdb14b05d03 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix comment
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 22 Jul 2022 16:31:18 +0200
parents af9e52e97154
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1268
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 classdef Elastic2dVariableAnisotropicMixedStencil < scheme.Scheme
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 % Discretizes the elastic wave equation:
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5 % opSet should be cell array of opSets, one per dimension. This
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 % is useful if we have periodic BC in one direction.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 % Assumes fully compatible operators
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9 properties
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 m % Number of points in each direction, possibly a vector
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 h % Grid spacing
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13 grid
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14 dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16 order % Order of accuracy for the approximation
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18 % Diagonal matrices for variable coefficients
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 RHO, RHOi, RHOi_kron % Density
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 C, C_D1, C_D2 % Elastic stiffness tensor, C = C_D1 + C_D2.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 D % Total operator
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 D1 % First derivatives
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24 % D2 % Second derivatives
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 % Boundary operators in cell format, used for BC
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 T_w, T_e, T_s, T_n
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29 % Traction operators
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 tau_w, tau_e, tau_s, tau_n % Return vector field
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 % Inner products
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35 H, Hi, Hi_kron, H_1D
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 % Boundary inner products (for scalar field)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 H_w, H_e, H_s, H_n
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 % Boundary restriction operators
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 % E{i}^T picks out component i
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 E
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 h11 % First entry in norm matrix
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 methods
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56 % Uses D1*D1 for the C_D1 part of the stiffness tensor C
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 % Uses narrow D2 whenever possible for the C_D2 part of C
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 % The coefficients can either be function handles or grid functions
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 function obj = Elastic2dVariableAnisotropicMixedStencil(g, order, rho, C_D1, C_D2, opSet)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 default_arg('rho', @(x,y) 0*x+1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 default_arg('optFlag', false);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 dim = 2;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65 C_D1_default = cell(dim,dim,dim,dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 C_D2_default = cell(dim,dim,dim,dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 for i = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 for j = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 for k = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70 for l = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 C_D1_default{i,j,k,l} = @(x,y) 0*x + 0;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 C_D2_default{i,j,k,l} = @(x,y) 0*x + 1;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 default_arg('C_D1', C_D1_default);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78 default_arg('C_D2', C_D2_default);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 assert(isa(g, 'grid.Cartesian'))
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 if isa(rho, 'function_handle')
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 rho = grid.evalOn(g, rho);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 C_mat = cell(dim,dim,dim,dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 C_D1_mat = cell(dim,dim,dim,dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 C_D2_mat = cell(dim,dim,dim,dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88 for i = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 for j = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 for k = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 for l = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 if isa(C_D1{i,j,k,l}, 'function_handle')
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 C_D1{i,j,k,l} = grid.evalOn(g, C_D1{i,j,k,l});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95 if isa(C_D2{i,j,k,l}, 'function_handle')
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 C_D2{i,j,k,l} = grid.evalOn(g, C_D2{i,j,k,l});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 C_D1_mat{i,j,k,l} = spdiag(C_D1{i,j,k,l});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 C_D2_mat{i,j,k,l} = spdiag(C_D2{i,j,k,l});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 C_mat{i,j,k,l} = C_D1_mat{i,j,k,l} + C_D2_mat{i,j,k,l};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 obj.C = C_mat;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 obj.C_D1 = C_D1_mat;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 obj.C_D2 = C_D2_mat;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 m = g.size();
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 m_tot = g.N();
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 lim = g.lim;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
112 if isempty(lim)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
113 x = g.x;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114 lim = cell(length(x),1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 for i = 1:length(x)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 lim{i} = {min(x{i}), max(x{i})};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 % 1D operators
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 ops = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 h = zeros(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 for i = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 ops{i} = opSet{i}(m(i), lim{i}, order);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 h(i) = ops{i}.h;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128 % Borrowing constants
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129 for i = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 obj.h11{i} = h(i)*ops{i}.borrowing.H11;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 I = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134 D1 = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135 D2 = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136 H = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137 Hi = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 e_0 = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139 e_m = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140 d1_0 = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141 d1_m = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143 for i = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144 I{i} = speye(m(i));
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145 D1{i} = ops{i}.D1;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 D2{i} = ops{i}.D2;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 H{i} = ops{i}.H;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148 Hi{i} = ops{i}.HI;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149 e_0{i} = ops{i}.e_l;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
150 e_m{i} = ops{i}.e_r;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
151 d1_0{i} = ops{i}.d1_l;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
152 d1_m{i} = ops{i}.d1_r;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155 %====== Assemble full operators ========
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
156 I_dim = speye(dim, dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
157 RHO = spdiag(rho);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
158 obj.RHO = RHO;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 obj.RHOi = inv(RHO);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
160 obj.RHOi_kron = kron(obj.RHOi, I_dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
161
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
162 obj.D1 = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
163 D2_temp = cell(dim,dim,dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
164
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
165 % D1
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
166 obj.D1{1} = kron(D1{1},I{2});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
167 obj.D1{2} = kron(I{1},D1{2});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
169 % Boundary restriction operators
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
170 e_l = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
171 e_r = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
172 e_l{1} = kron(e_0{1}, I{2});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173 e_l{2} = kron(I{1}, e_0{2});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174 e_r{1} = kron(e_m{1}, I{2});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 e_r{2} = kron(I{1}, e_m{2});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
176
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
177 e_scalar_w = e_l{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178 e_scalar_e = e_r{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179 e_scalar_s = e_l{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 e_scalar_n = e_r{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
182 e_w = kron(e_scalar_w, I_dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
183 e_e = kron(e_scalar_e, I_dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184 e_s = kron(e_scalar_s, I_dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
185 e_n = kron(e_scalar_n, I_dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187 % E{i}^T picks out component i.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
188 E = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
189 I = speye(m_tot,m_tot);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 for i = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191 e = sparse(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192 e(i) = 1;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
193 E{i} = kron(I,e);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
194 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195 obj.E = E;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
197 e1_w = (e_scalar_w'*E{1}')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198 e1_e = (e_scalar_e'*E{1}')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 e1_s = (e_scalar_s'*E{1}')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200 e1_n = (e_scalar_n'*E{1}')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 e2_w = (e_scalar_w'*E{2}')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 e2_e = (e_scalar_e'*E{2}')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 e2_s = (e_scalar_s'*E{2}')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 e2_n = (e_scalar_n'*E{2}')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 % D2
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 switch order
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210 case 2
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 width = 3;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212 case 4
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
213 width = 5;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214 case 6
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 width = 7;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 for j = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
218 for k = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219 for l = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 ind = grid.funcToMatrix(g, 1:m_tot);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226 k = 1;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
227 for r = 1:m(2)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228 p = ind(:,r);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 for j = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 for l = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231 coeff = C_D2{k,j,k,l};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
232 D_kk = D2{1}(coeff(p));
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
233 D2_temp{j,k,l}(p,p) = D_kk;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
234 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
235 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
236 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
237
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
238 k = 2;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
239 for r = 1:m(1)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 p = ind(r,:);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241 for j = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
242 for l = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
243 coeff = C_D2{k,j,k,l};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
244 D_kk = D2{2}(coeff(p));
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
245 D2_temp{j,k,l}(p,p) = D_kk;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
246 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
247 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
249
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
250 % Quadratures
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
251 obj.H = kron(H{1},H{2});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
252 obj.Hi = inv(obj.H);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253 obj.H_w = H{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
254 obj.H_e = H{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 obj.H_s = H{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256 obj.H_n = H{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
257 obj.H_1D = {H{1}, H{2}};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
258
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
259 % Differentiation matrix D (without SAT)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
260 D1 = obj.D1;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
261 D = sparse(dim*m_tot,dim*m_tot);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
262 for i = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
263 for j = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
264 for k = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
265 for l = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
266 if i == k
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
267 D = D + E{j}*D2_temp{j,k,l}*E{l}';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
268 D2_temp{j,k,l} = [];
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
269 else
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270 D = D + E{j}*D1{i}*C_D2_mat{i,j,k,l}*D1{k}*E{l}';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
271 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
272 D = D + E{j}*D1{i}*C_D1_mat{i,j,k,l}*D1{k}*E{l}';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
273 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
274 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
275 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
276 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
277 clear D2_temp;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
278 D = obj.RHOi_kron*D;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
279 obj.D = D;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
280 clear D;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
281 %=========================================%'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
282
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
283 % Numerical traction operators for BC.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
284 %
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
285 % Formula at boundary j: % tau^{j}_i = sum_l T^{j}_{il} u_l
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
286 %
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
287 T_l = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
288 T_r = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
289 tau_l = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
290 tau_r = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
291
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
292 D1 = obj.D1;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
293
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
294 % Boundary j
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
295 for j = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
296 T_l{j} = cell(dim,dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
297 T_r{j} = cell(dim,dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
298 tau_l{j} = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
299 tau_r{j} = cell(dim,1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
300
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
301 [~, n_l] = size(e_l{j});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
302 [~, n_r] = size(e_r{j});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
303
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
304 % Traction component i
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
305 for i = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
306 tau_l{j}{i} = sparse(dim*m_tot, n_l);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
307 tau_r{j}{i} = sparse(dim*m_tot, n_r);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
308
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
309 % Displacement component l
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
310 for l = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
311 T_l{j}{i,l} = sparse(m_tot, n_l);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
312 T_r{j}{i,l} = sparse(m_tot, n_r);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
313
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
314 % Derivative direction k
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
315 for k = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
316 T_l{j}{i,l} = T_l{j}{i,l} ...
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
317 - (e_l{j}'*C_mat{j,i,k,l}*D1{k})';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
318 T_r{j}{i,l} = T_r{j}{i,l} ...
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
319 + (e_r{j}'*C_mat{j,i,k,l}*D1{k})';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
320 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
321 tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
322 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
323 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
324 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
325 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
326
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
327 % Traction tensors, T_ij
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
328 obj.T_w = T_l{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
329 obj.T_e = T_r{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
330 obj.T_s = T_l{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
331 obj.T_n = T_r{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
332
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
333 % Restriction operators
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
334 obj.e_w = e_w;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
335 obj.e_e = e_e;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
336 obj.e_s = e_s;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
337 obj.e_n = e_n;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
338
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
339 obj.e1_w = e1_w;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
340 obj.e1_e = e1_e;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
341 obj.e1_s = e1_s;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
342 obj.e1_n = e1_n;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
343
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
344 obj.e2_w = e2_w;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
345 obj.e2_e = e2_e;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
346 obj.e2_s = e2_s;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
347 obj.e2_n = e2_n;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
348
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
349 obj.e_scalar_w = e_scalar_w;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
350 obj.e_scalar_e = e_scalar_e;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 obj.e_scalar_s = e_scalar_s;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352 obj.e_scalar_n = e_scalar_n;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
353
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
354 % First component of traction
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
355 obj.tau1_w = tau_l{1}{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
356 obj.tau1_e = tau_r{1}{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
357 obj.tau1_s = tau_l{2}{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
358 obj.tau1_n = tau_r{2}{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
359
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
360 % Second component of traction
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
361 obj.tau2_w = tau_l{1}{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
362 obj.tau2_e = tau_r{1}{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
363 obj.tau2_s = tau_l{2}{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
364 obj.tau2_n = tau_r{2}{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
365
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
366 % Traction vectors
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
367 obj.tau_w = (e_w'*e1_w*obj.tau1_w')' + (e_w'*e2_w*obj.tau2_w')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
368 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
369 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
370 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
371
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
372 % Kroneckered norms and coefficients
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
373 obj.Hi_kron = kron(obj.Hi, I_dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
374
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
375 % Misc.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
376 obj.m = m;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377 obj.h = h;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
378 obj.order = order;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
379 obj.grid = g;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
380 obj.dim = dim;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
381
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
382 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
383
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
384
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
385 % Closure functions return the operators applied to the own domain to close the boundary
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
386 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
387 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
388 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
389 % on the first component. Can also be e.g.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
390 % {'normal', 'd'} or {'tangential', 't'} for conditions on
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
391 % tangential/normal component.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
392 % data is a function returning the data that should be applied at the boundary.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
393 % neighbour_scheme is an instance of Scheme that should be interfaced to.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
394 % neighbour_boundary is a string specifying which boundary to interface to.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
395
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
396 % For displacement bc:
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
397 % bc = {comp, 'd', dComps},
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
398 % where
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
399 % dComps = vector of components with displacement BC. Default: 1:dim.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
400 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
401 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
402 default_arg('tuning', 1.0);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
403
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
404 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
405 comp = bc{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
406 type = bc{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
407 if ischar(comp)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
408 comp = obj.getComponent(comp, boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
409 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
410
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
411 e = obj.getBoundaryOperatorForScalarField('e', boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
412 tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
413 T = obj.getBoundaryTractionOperator(boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
414 h11 = obj.getBorrowing(boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
415 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
416 nu = obj.getNormal(boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
417
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
418 E = obj.E;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
419 Hi = obj.Hi;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
420 RHOi = obj.RHOi;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
421 C = obj.C;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
422
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
423 dim = obj.dim;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
424 m_tot = obj.grid.N();
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
425
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
426 % Preallocate
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
427 [~, col] = size(tau);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
428 closure = sparse(dim*m_tot, dim*m_tot);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
429 penalty = sparse(dim*m_tot, col);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
430
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
431 j = comp;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
432 switch type
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
433
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
434 % Dirichlet boundary condition
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
435 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
436
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
437 if numel(bc) >= 3
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
438 dComps = bc{3};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
439 else
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
440 dComps = 1:dim;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
441 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
442
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
443 % Loops over components that Dirichlet penalties end up on
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
444 % Y: symmetrizing part of penalty
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
445 % Z: symmetric part of penalty
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
446 % X = Y + Z.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
447
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
448 % Nonsymmetric part goes on all components to
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
449 % yield traction in discrete energy rate
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
450 for i = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
451 Y = T{j,i}';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
452 X = e*Y;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
453 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
454 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
455 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
456
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
457 % Symmetric part only required on components with displacement BC.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
458 % (Otherwise it's not symmetric.)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
459 for i = dComps
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
460 Z = sparse(m_tot, m_tot);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
461 for l = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
462 for k = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
463 Z = Z + nu(l)*C{l,i,k,j}*nu(k);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
464 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
465 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
466 Z = -tuning*dim/h11*Z;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
467 X = Z;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
468 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
469 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
470 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
471
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
472 % Free boundary condition
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
473 case {'F','f','Free','free','traction','Traction','t','T'}
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
474 closure = closure - E{j}*RHOi*Hi*e*H_gamma*tau';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
475 penalty = penalty + E{j}*RHOi*Hi*e*H_gamma;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
476
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
477 % Unknown boundary condition
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
478 otherwise
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
479 error('No such boundary condition: type = %s',type);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
480 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
481 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
482
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
483 % type Struct that specifies the interface coupling.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
484 % Fields:
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
485 % -- tuning: penalty strength, defaults to 1.0
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
486 % -- interpolation: type of interpolation, default 'none'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
487 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
488
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
489 defaultType.tuning = 1.0;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
490 defaultType.interpolation = 'none';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
491 default_struct('type', defaultType);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
492
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
493 switch type.interpolation
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
494 case {'none', ''}
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
495 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
496 case {'op','OP'}
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
497 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
498 otherwise
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
499 error('Unknown type of interpolation: %s ', type.interpolation);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
500 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
501 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
502
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
503 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
504 tuning = type.tuning;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
505
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
506 % u denotes the solution in the own domain
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
507 % v denotes the solution in the neighbour domain
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
508
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
509 u = obj;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
510 v = neighbour_scheme;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
511
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
512 % Operators, u side
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
513 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
514 tau_u = u.getBoundaryOperator('tau', boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
515 h11_u = u.getBorrowing(boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
516 nu_u = u.getNormal(boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
517
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
518 E_u = u.E;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
519 C_u = u.C;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
520 m_tot_u = u.grid.N();
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
521
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
522 % Operators, v side
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
523 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
524 tau_v = v.getBoundaryOperator('tau', neighbour_boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
525 h11_v = v.getBorrowing(neighbour_boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
526 nu_v = v.getNormal(neighbour_boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
527
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
528 E_v = v.E;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
529 C_v = v.C;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
530 m_tot_v = v.grid.N();
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
531
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
532 % Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
533 flipFlag = false;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
534 e_v_flip = e_v;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
535 if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) || ...
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
536 (strcmp(boundary,'e') && strcmp(neighbour_boundary,'s')) || ...
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
537 (strcmp(boundary,'w') && strcmp(neighbour_boundary,'n')) || ...
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
538 (strcmp(boundary,'n') && strcmp(neighbour_boundary,'w')) || ...
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
539 (strcmp(boundary,'s') && strcmp(neighbour_boundary,'s')) || ...
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
540 (strcmp(boundary,'n') && strcmp(neighbour_boundary,'n')) || ...
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
541 (strcmp(boundary,'w') && strcmp(neighbour_boundary,'w')) || ...
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
542 (strcmp(boundary,'e') && strcmp(neighbour_boundary,'e'))
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
543
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
544 flipFlag = true;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
545 e_v_flip = fliplr(e_v);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
546
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
547 t1 = tau_v(:,1:2:end-1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
548 t2 = tau_v(:,2:2:end);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
549
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
550 t1 = fliplr(t1);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
551 t2 = fliplr(t2);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
552
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
553 tau_v(:,1:2:end-1) = t1;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
554 tau_v(:,2:2:end) = t2;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
555 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
556
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
557 % Operators that are only required for own domain
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
558 Hi = u.Hi_kron;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
559 RHOi = u.RHOi_kron;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
560 e_kron = u.getBoundaryOperator('e', boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
561 T_u = u.getBoundaryTractionOperator(boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
562
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
563 % Shared operators
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
564 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
565 H_gamma_kron = u.getBoundaryQuadrature(boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
566 dim = u.dim;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
567
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
568 % Preallocate
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
569 [~, m_int] = size(H_gamma);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
570 closure = sparse(dim*m_tot_u, dim*m_tot_u);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
571 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
572
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
573 % ---- Continuity of displacement ------
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
574
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
575 % Y: symmetrizing part of penalty
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
576 % Z: symmetric part of penalty
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
577 % X = Y + Z.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
578
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
579 % Loop over components to couple across interface
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
580 for j = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
581
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
582 % Loop over components that penalties end up on
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
583 for i = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
584 Y = 1/2*T_u{j,i}';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
585 Z_u = sparse(m_int, m_int);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
586 Z_v = sparse(m_int, m_int);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
587 for l = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
588 for k = 1:dim
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
589 Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
590 Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
591 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
592 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
593
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
594 if flipFlag
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
595 Z_v = rot90(Z_v,2);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
596 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
597
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
598 Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
599 X = Y + Z*e_u';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
600 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
601 penalty = penalty - E_u{i}*X'*H_gamma*e_v_flip'*E_v{j}';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
602
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
603 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
604 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
605
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
606 % ---- Continuity of traction ------
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
607 closure = closure - 1/2*e_kron*H_gamma_kron*tau_u';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
608 penalty = penalty - 1/2*e_kron*H_gamma_kron*tau_v';
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
609
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
610 % ---- Multiply by inverse of density x quadraure ----
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
611 closure = RHOi*Hi*closure;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
612 penalty = RHOi*Hi*penalty;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
613
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
614 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
615
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
616 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
617 error('Non-conforming interfaces not implemented yet.');
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
618 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
619
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
620 % Returns the component number that is the tangential/normal component
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
621 % at the specified boundary
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
622 function comp = getComponent(obj, comp_str, boundary)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
623 assertIsMember(comp_str, {'normal', 'tangential'});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
624 assertIsMember(boundary, {'w', 'e', 's', 'n'});
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
625
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
626 switch boundary
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
627 case {'w', 'e'}
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
628 switch comp_str
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
629 case 'normal'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
630 comp = 1;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
631 case 'tangential'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
632 comp = 2;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
633 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
634 case {'s', 'n'}
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
635 switch comp_str
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
636 case 'normal'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
637 comp = 2;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
638 case 'tangential'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
639 comp = 1;
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
640 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
641 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
642 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
643
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
644 % Returns h11 for the boundary specified by the string boundary.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
645 % op -- string
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
646 function h11 = getBorrowing(obj, boundary)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
647 assertIsMember(boundary, {'w', 'e', 's', 'n'})
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
648
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
649 switch boundary
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
650 case {'w','e'}
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
651 h11 = obj.h11{1};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
652 case {'s', 'n'}
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
653 h11 = obj.h11{2};
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
654 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
655 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
656
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
657 % Returns the outward unit normal vector for the boundary specified by the string boundary.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
658 function nu = getNormal(obj, boundary)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
659 assertIsMember(boundary, {'w', 'e', 's', 'n'})
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
660
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
661 switch boundary
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
662 case 'w'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
663 nu = [-1,0];
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
664 case 'e'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
665 nu = [1,0];
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
666 case 's'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
667 nu = [0,-1];
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
668 case 'n'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
669 nu = [0,1];
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
670 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
671 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
672
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
673 % Returns the boundary operator op for the boundary specified by the string boundary.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
674 % op -- string
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
675 function o = getBoundaryOperator(obj, op, boundary)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
676 assertIsMember(boundary, {'w', 'e', 's', 'n'})
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
677 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'})
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
678
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
679 switch op
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
680 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
681 o = obj.([op, '_', boundary]);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
682 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
683
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
684 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
685
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
686 % Returns the boundary operator op for the boundary specified by the string boundary.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
687 % op -- string
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
688 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
689 assertIsMember(boundary, {'w', 'e', 's', 'n'})
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
690 assertIsMember(op, {'e'})
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
691
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
692 switch op
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
693
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
694 case 'e'
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
695 o = obj.(['e_scalar', '_', boundary]);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
696 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
697
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
698 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
699
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
700 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
701 % Formula: tau_i = T_ij u_j
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
702 % op -- string
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
703 function T = getBoundaryTractionOperator(obj, boundary)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
704 assertIsMember(boundary, {'w', 'e', 's', 'n'})
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
705
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
706 T = obj.(['T', '_', boundary]);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
707 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
708
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
709 % Returns square boundary quadrature matrix, of dimension
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
710 % corresponding to the number of boundary unknowns
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
711 %
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
712 % boundary -- string
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
713 function H = getBoundaryQuadrature(obj, boundary)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
714 assertIsMember(boundary, {'w', 'e', 's', 'n'})
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
715
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
716 H = obj.getBoundaryQuadratureForScalarField(boundary);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
717 I_dim = speye(obj.dim, obj.dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
718 H = kron(H, I_dim);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
719 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
720
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
721 % Returns square boundary quadrature matrix, of dimension
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
722 % corresponding to the number of boundary grid points
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
723 %
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
724 % boundary -- string
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
725 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
726 assertIsMember(boundary, {'w', 'e', 's', 'n'})
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
727
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
728 H_b = obj.(['H_', boundary]);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
729 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
730
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
731 function N = size(obj)
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
732 N = obj.dim*prod(obj.m);
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
733 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
734 end
af9e52e97154 Add Elastic2dVariableAnisotropicMixedStencil
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
735 end