Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropicUpwind.m @ 1305:b5907140c069 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up hollow D2 generation by using precomputed coefficients stored in mat files
author Martin Almquist <malmquist@stanford.edu>
date Mon, 13 Jul 2020 14:21:41 -0700
parents 8ff3a95ad7cc
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1227
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 classdef Elastic2dCurvilinearAnisotropicUpwind < 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 % in curvilinear coordinates.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 % 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
7 % 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
8 % Assumes fully compatible operators.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 properties
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 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
12 h % Grid spacing
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14 grid
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 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
18
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 % Diagonal matrices for variable coefficients
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 J, Ji
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21 RHO % Density
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 C % Elastic stiffness tensor
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24 D % Total operator
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 Dx, Dy % Physical derivatives
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 sigma % Cell matrix of physical stress operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28 n_w, n_e, n_s, n_n % Physical normals
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 % 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
31 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
32
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 % Traction operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 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
35 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
36 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
37
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 % Inner products
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39 H
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 % Boundary inner products (for scalar field)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 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
43
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 % Surface Jacobian vectors
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45 s_w, s_e, s_s, s_n
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 % Boundary restriction operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 en_w, en_e, en_s, en_n % Act on vector field, return normal component
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 % 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
55 E
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 % Elastic2dVariableAnisotropic object for reference domain
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 refObj
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61 methods
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 % 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
64 % 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
65 function obj = Elastic2dCurvilinearAnisotropicUpwind(g, order, rho, C, opSet, optFlag)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 default_arg('rho', @(x,y) 0*x+1);
1270
9a7cac4202d0 In CruvilinearAnisotropicUpwind: Make D1Upwind default opSet. Use D2Standard for metric derivatives. Bugfix global quadrature.
Martin Almquist <malmquist@stanford.edu>
parents: 1227
diff changeset
67 default_arg('opSet',{@sbp.D1Upwind, @sbp.D1Upwind});
1227
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 default_arg('optFlag', false);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 dim = 2;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 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
72 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 for l = 1:dim
1272
15865fbda16e Set default stiffness tensors to 0 in CurvilinearAnisotropic schemes
Martin Almquist <malmquist@stanford.edu>
parents: 1270
diff changeset
76 C_default{i,j,k,l} = @(x,y) 0*x;
1227
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 end
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 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81 default_arg('C', C_default);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 assert(isa(g, 'grid.Curvilinear'));
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 if isa(rho, 'function_handle')
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 rho = grid.evalOn(g, rho);
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
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 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
90 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 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
95 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
96 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 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
98 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
101 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 obj.C = C_mat;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 m = g.size();
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 m_tot = g.N();
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 % 1D operators
1270
9a7cac4202d0 In CruvilinearAnisotropicUpwind: Make D1Upwind default opSet. Use D2Standard for metric derivatives. Bugfix global quadrature.
Martin Almquist <malmquist@stanford.edu>
parents: 1227
diff changeset
108 opSetMetric = {@sbp.D2VariableCompatible, @sbp.D2VariableCompatible};
1275
8ff3a95ad7cc Change metric computation in CurvilinearUpwind so that one can use odd orders for the diffOp
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
109 orderMetric = ceil(order/2)*2;
1227
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 m_u = m(1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
111 m_v = m(2);
1275
8ff3a95ad7cc Change metric computation in CurvilinearUpwind so that one can use odd orders for the diffOp
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
112 ops_u = opSetMetric{1}(m_u, {0, 1}, orderMetric);
8ff3a95ad7cc Change metric computation in CurvilinearUpwind so that one can use odd orders for the diffOp
Martin Almquist <malmquist@stanford.edu>
parents: 1272
diff changeset
113 ops_v = opSetMetric{2}(m_v, {0, 1}, orderMetric);
1227
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
114
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 h_u = ops_u.h;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116 h_v = ops_v.h;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118 I_u = speye(m_u);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119 I_v = speye(m_v);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 D1_u = ops_u.D1;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 H_u = ops_u.H;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 Hi_u = ops_u.HI;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
124 e_l_u = ops_u.e_l;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
125 e_r_u = ops_u.e_r;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
126 d1_l_u = ops_u.d1_l;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
127 d1_r_u = ops_u.d1_r;
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 D1_v = ops_v.D1;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 H_v = ops_v.H;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 Hi_v = ops_v.HI;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
132 e_l_v = ops_v.e_l;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
133 e_r_v = ops_v.e_r;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
134 d1_l_v = ops_v.d1_l;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
135 d1_r_v = ops_v.d1_r;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
136
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
137
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 % Logical operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139 Du = kr(D1_u,I_v);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
140 Dv = kr(I_u,D1_v);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
141
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
142 e_w = kr(e_l_u,I_v);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
143 e_e = kr(e_r_u,I_v);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
144 e_s = kr(I_u,e_l_v);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
145 e_n = kr(I_u,e_r_v);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 % Metric coefficients
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148 coords = g.points();
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149 x = coords(:,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
150 y = coords(:,2);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
151
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
152 x_u = Du*x;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 x_v = Dv*x;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154 y_u = Du*y;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155 y_v = Dv*y;
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 J = x_u.*y_v - x_v.*y_u;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
158
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 K = cell(dim, dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
160 K{1,1} = y_v./J;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
161 K{1,2} = -y_u./J;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
162 K{2,1} = -x_v./J;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
163 K{2,2} = x_u./J;
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 % Physical derivatives
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
166 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
167 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
169 % Wrap around Aniosotropic Cartesian
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
170 rho_tilde = J.*rho;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
171
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
172 PHI = cell(dim,dim,dim,dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
176 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
177 PHI{i,j,k,l} = 0*C{i,j,k,l};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178 for m = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179 for n = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
180 PHI{i,j,k,l} = PHI{i,j,k,l} + J.*K{m,i}.*C{m,j,n,l}.*K{n,k};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
181 end
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 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
184 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
185 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
186 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
188 gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1});
1270
9a7cac4202d0 In CruvilinearAnisotropicUpwind: Make D1Upwind default opSet. Use D2Standard for metric derivatives. Bugfix global quadrature.
Martin Almquist <malmquist@stanford.edu>
parents: 1227
diff changeset
189 refObj = scheme.Elastic2dVariableAnisotropicUpwind(gRef, order, rho_tilde, PHI, opSet);
1227
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191 %---- Set object properties ------
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192 obj.RHO = spdiag(rho);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
193
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
194 % Volume quadrature
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195 obj.J = spdiag(J);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
196 obj.Ji = spdiag(1./J);
1270
9a7cac4202d0 In CruvilinearAnisotropicUpwind: Make D1Upwind default opSet. Use D2Standard for metric derivatives. Bugfix global quadrature.
Martin Almquist <malmquist@stanford.edu>
parents: 1227
diff changeset
197 obj.H = obj.J*refObj.H;
1227
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 % Boundary quadratures
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
200 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
201 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
202 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
203 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 obj.s_w = s_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 obj.s_e = s_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206 obj.s_s = s_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207 obj.s_n = s_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
209 obj.H_w = H_v*spdiag(s_w);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
210 obj.H_e = H_v*spdiag(s_e);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
211 obj.H_s = H_u*spdiag(s_s);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
212 obj.H_n = H_u*spdiag(s_n);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
213
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214 % Restriction operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 obj.e_w = refObj.e_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 obj.e_e = refObj.e_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 obj.e_s = refObj.e_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
218 obj.e_n = refObj.e_n;
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 % Adapt things from reference object
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221 obj.D = refObj.D;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
222 obj.E = refObj.E;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
223
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 obj.e1_w = refObj.e1_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 obj.e1_e = refObj.e1_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226 obj.e1_s = refObj.e1_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
227 obj.e1_n = refObj.e1_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 obj.e2_w = refObj.e2_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 obj.e2_e = refObj.e2_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231 obj.e2_s = refObj.e2_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
232 obj.e2_n = refObj.e2_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
233
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
234 obj.e_scalar_w = refObj.e_scalar_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
235 obj.e_scalar_e = refObj.e_scalar_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
236 obj.e_scalar_s = refObj.e_scalar_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
237 obj.e_scalar_n = refObj.e_scalar_n;
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 e1_w = obj.e1_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 e1_e = obj.e1_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241 e1_s = obj.e1_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
242 e1_n = obj.e1_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
243
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
244 e2_w = obj.e2_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
245 e2_e = obj.e2_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
246 e2_s = obj.e2_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
247 e2_n = obj.e2_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
249 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
250 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
251 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
252 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
254 obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
255 obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
256 obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
257 obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
258
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
259 obj.tau_w = (refObj.e_w'*obj.e1_w*obj.tau1_w')' + (refObj.e_w'*obj.e2_w*obj.tau2_w')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
260 obj.tau_e = (refObj.e_e'*obj.e1_e*obj.tau1_e')' + (refObj.e_e'*obj.e2_e*obj.tau2_e')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
261 obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
262 obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
263
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
264 % Physical normals
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
265 e_w = obj.e_scalar_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
266 e_e = obj.e_scalar_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
267 e_s = obj.e_scalar_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
268 e_n = obj.e_scalar_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
269
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270 e_w_vec = obj.e_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
271 e_e_vec = obj.e_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
272 e_s_vec = obj.e_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
273 e_n_vec = obj.e_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
274
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
275 nu_w = [-1,0];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
276 nu_e = [1,0];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
277 nu_s = [0,-1];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
278 nu_n = [0,1];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
279
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
280 obj.n_w = cell(2,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
281 obj.n_e = cell(2,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
282 obj.n_s = cell(2,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
283 obj.n_n = cell(2,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
284
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
285 n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2)));
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
286 n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2)));
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
287 obj.n_w{1} = spdiag(n_w_1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
288 obj.n_w{2} = spdiag(n_w_2);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
289
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
290 n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2)));
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
291 n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2)));
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
292 obj.n_e{1} = spdiag(n_e_1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
293 obj.n_e{2} = spdiag(n_e_2);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
294
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
295 n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2)));
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
296 n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2)));
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
297 obj.n_s{1} = spdiag(n_s_1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
298 obj.n_s{2} = spdiag(n_s_2);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
299
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
300 n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2)));
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
301 n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2)));
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
302 obj.n_n{1} = spdiag(n_n_1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
303 obj.n_n{2} = spdiag(n_n_2);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
304
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
305 % Operators that extract the normal component
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
306 obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
307 obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
308 obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
309 obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
310
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
311 % Stress operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
312 sigma = cell(dim, dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
313 D1 = {obj.Dx, obj.Dy};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
314 E = obj.E;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
315 N = length(obj.RHO);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
316 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
317 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
318 sigma{i,j} = sparse(N,2*N);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
319 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
320 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
321 sigma{i,j} = sigma{i,j} + obj.C{i,j,k,l}*D1{k}*E{l}';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
322 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
323 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
324 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
325 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
326 obj.sigma = sigma;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
327
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
328 % Misc.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
329 obj.refObj = refObj;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
330 obj.m = refObj.m;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
331 obj.h = refObj.h;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
332 obj.order = order;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
333 obj.grid = g;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
334 obj.dim = dim;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
335
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
336 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
337
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
338
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
339 % 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
340 % 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
341 % 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
342 % 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
343 % 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
344 % {'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
345 % tangential/normal component.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
346 % 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
347 % 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
348 % 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
349
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
350 % For displacement bc:
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 % bc = {comp, 'd', dComps},
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352 % where
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
353 % 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
354 % 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
355 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
356 default_arg('tuning', 1.0);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
357 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
358
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
359 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
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 type = bc{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
362
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
363 switch type
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
364 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
365 s = obj.(['s_' boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
366 s = spdiag(s);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
367 penalty = penalty*s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
368 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
369 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
370
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
371 % 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
372 % Fields:
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
373 % -- 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
374 % -- interpolation: type of interpolation, default 'none'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
375 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
376
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377 defaultType.tuning = 1.0;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
378 defaultType.interpolation = 'none';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
379 default_struct('type', defaultType);
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 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
382 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
383
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
384 % 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
385 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
386 function h11 = getBorrowing(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
387 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
388
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
389 switch boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
390 case {'w','e'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
391 h11 = obj.refObj.h11{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
392 case {'s', 'n'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
393 h11 = obj.refObj.h11{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
394 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
395 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
396
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
397 % 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
398 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
399 function n = getNormal(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
400 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
401
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
402 n = obj.(['n_' boundary]);
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
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
405 % 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
406 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
407 function o = getBoundaryOperator(obj, op, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
408 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
409 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
410
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
411 o = obj.([op, '_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
412
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
413 end
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 % 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
416 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
417 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
418 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
419 assertIsMember(op, {'e'})
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 switch op
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
422
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
423 case 'e'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
424 o = obj.(['e_scalar', '_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
425 end
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 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
428
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
429 % 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
430 % 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
431 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
432 function T = getBoundaryTractionOperator(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
433 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
434
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
435 T = obj.(['T', '_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
436 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
437
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
438 % 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
439 % 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
440 %
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
441 % boundary -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
442 function H = getBoundaryQuadrature(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
443 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
444
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
445 H = obj.getBoundaryQuadratureForScalarField(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
446 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
447 H = kron(H, I_dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
448 end
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 % 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
451 % 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
452 %
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
453 % boundary -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
454 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
455 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
456
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
457 H_b = obj.(['H_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
458 end
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 function N = size(obj)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
461 N = obj.dim*prod(obj.m);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
462 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
463 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
464 end