Mercurial > repos > public > sbplib
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 |
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 |