annotate +scheme/Elastic2dVariableAnisotropicUpwind.m @ 1253:89dad61cad22 feature/poroelastic

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