Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropicUpwind.m @ 1272:15865fbda16e feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Set default stiffness tensors to 0 in CurvilinearAnisotropic schemes
author Martin Almquist <malmquist@stanford.edu>
date Sun, 31 May 2020 21:59:09 -0700
parents 9a7cac4202d0
children 8ff3a95ad7cc
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};
1227
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
109 m_u = m(1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
110 m_v = m(2);
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
111 ops_u = opSetMetric{1}(m_u, {0, 1}, order);
9a7cac4202d0 In CruvilinearAnisotropicUpwind: Make D1Upwind default opSet. Use D2Standard for metric derivatives. Bugfix global quadrature.
Martin Almquist <malmquist@stanford.edu>
parents: 1227
diff changeset
112 ops_v = opSetMetric{2}(m_v, {0, 1}, order);
1227
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 h_u = ops_u.h;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
115 h_v = ops_v.h;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
116
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
117 I_u = speye(m_u);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
118 I_v = speye(m_v);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
119
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
120 D1_u = ops_u.D1;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
121 H_u = ops_u.H;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
122 Hi_u = ops_u.HI;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
123 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
124 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
125 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
126 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
127
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
128 D1_v = ops_v.D1;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
129 H_v = ops_v.H;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
130 Hi_v = ops_v.HI;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
131 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
132 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
133 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
134 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
135
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 % Logical operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
138 Du = kr(D1_u,I_v);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
139 Dv = kr(I_u,D1_v);
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 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
142 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
143 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
144 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
145
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
146 % Metric coefficients
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
147 coords = g.points();
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
148 x = coords(:,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
149 y = coords(:,2);
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 x_u = Du*x;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
152 x_v = Dv*x;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
153 y_u = Du*y;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
154 y_v = Dv*y;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
155
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
156 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
157
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
158 K = cell(dim, dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
159 K{1,1} = y_v./J;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
160 K{1,2} = -y_u./J;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
161 K{2,1} = -x_v./J;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
162 K{2,2} = x_u./J;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
163
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
164 % Physical derivatives
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
165 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
166 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
167
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
168 % Wrap around Aniosotropic Cartesian
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
169 rho_tilde = J.*rho;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
170
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
171 PHI = cell(dim,dim,dim,dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
172 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
173 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
174 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
175 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
176 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
177 for m = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
178 for n = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
179 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
180 end
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
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
187 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
188 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
189
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
190 %---- Set object properties ------
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
191 obj.RHO = spdiag(rho);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
192
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
193 % Volume quadrature
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
194 obj.J = spdiag(J);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
195 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
196 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
197
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
198 % Boundary quadratures
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
199 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
200 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
201 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
202 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
203 obj.s_w = s_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
204 obj.s_e = s_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
205 obj.s_s = s_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
206 obj.s_n = s_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
207
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
208 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
209 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
210 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
211 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
212
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
213 % Restriction operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
214 obj.e_w = refObj.e_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
215 obj.e_e = refObj.e_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
216 obj.e_s = refObj.e_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
217 obj.e_n = refObj.e_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
218
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
219 % Adapt things from reference object
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
220 obj.D = refObj.D;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
221 obj.E = refObj.E;
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 obj.e1_w = refObj.e1_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
224 obj.e1_e = refObj.e1_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
225 obj.e1_s = refObj.e1_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
226 obj.e1_n = refObj.e1_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
227
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
228 obj.e2_w = refObj.e2_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
229 obj.e2_e = refObj.e2_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
230 obj.e2_s = refObj.e2_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
231 obj.e2_n = refObj.e2_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
232
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
233 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
234 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
235 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
236 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
237
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
238 e1_w = obj.e1_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
239 e1_e = obj.e1_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
240 e1_s = obj.e1_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
241 e1_n = obj.e1_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
242
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
243 e2_w = obj.e2_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
244 e2_e = obj.e2_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
245 e2_s = obj.e2_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
246 e2_n = obj.e2_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
247
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
248 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
249 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
250 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
251 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
252
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
253 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
254 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
255 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
256 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
257
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
258 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
259 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
260 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
261 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
262
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
263 % Physical normals
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
264 e_w = obj.e_scalar_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
265 e_e = obj.e_scalar_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
266 e_s = obj.e_scalar_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
267 e_n = obj.e_scalar_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
268
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
269 e_w_vec = obj.e_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
270 e_e_vec = obj.e_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
271 e_s_vec = obj.e_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
272 e_n_vec = obj.e_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
273
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
274 nu_w = [-1,0];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
275 nu_e = [1,0];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
276 nu_s = [0,-1];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
277 nu_n = [0,1];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
278
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
279 obj.n_w = cell(2,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
280 obj.n_e = cell(2,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
281 obj.n_s = cell(2,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
282 obj.n_n = cell(2,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
283
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
284 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
285 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
286 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
287 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
288
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
289 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
290 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
291 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
292 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
293
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
294 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
295 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
296 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
297 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
298
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
299 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
300 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
301 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
302 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
303
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
304 % Operators that extract the normal component
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
305 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
306 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
307 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
308 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
309
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
310 % Stress operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
311 sigma = cell(dim, dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
312 D1 = {obj.Dx, obj.Dy};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
313 E = obj.E;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
314 N = length(obj.RHO);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
315 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
316 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
317 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
318 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
319 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
320 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
321 end
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 obj.sigma = sigma;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
326
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
327 % Misc.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
328 obj.refObj = refObj;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
329 obj.m = refObj.m;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
330 obj.h = refObj.h;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
331 obj.order = order;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
332 obj.grid = g;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
333 obj.dim = dim;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
334
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
335 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
336
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 % 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
339 % 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
340 % 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
341 % 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
342 % 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
343 % {'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
344 % tangential/normal component.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
345 % 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
346 % 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
347 % 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
348
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
349 % For displacement bc:
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
350 % bc = {comp, 'd', dComps},
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
351 % where
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
352 % 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
353 % 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
354 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
355 default_arg('tuning', 1.0);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
356 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
357
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
358 [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
359
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
360 type = bc{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
361
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
362 switch type
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
363 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
364 s = obj.(['s_' boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
365 s = spdiag(s);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
366 penalty = penalty*s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
367 end
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
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
370 % 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
371 % Fields:
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
372 % -- 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
373 % -- interpolation: type of interpolation, default 'none'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
374 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
375
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
376 defaultType.tuning = 1.0;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
377 defaultType.interpolation = 'none';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
378 default_struct('type', defaultType);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
379
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
380 [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
381 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
382
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
383 % 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
384 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
385 function h11 = getBorrowing(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
386 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
387
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
388 switch boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
389 case {'w','e'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
390 h11 = obj.refObj.h11{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
391 case {'s', 'n'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
392 h11 = obj.refObj.h11{2};
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 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
395
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
396 % 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
397 % 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
398 function n = getNormal(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
399 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
400
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
401 n = obj.(['n_' boundary]);
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
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
404 % 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
405 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
406 function o = getBoundaryOperator(obj, op, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
407 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
408 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
409
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
410 o = obj.([op, '_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
411
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
412 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
413
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
414 % 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
415 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
416 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
417 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
418 assertIsMember(op, {'e'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
419
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
420 switch op
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
421
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
422 case 'e'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
423 o = obj.(['e_scalar', '_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
424 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
425
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
426 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
427
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
428 % 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
429 % 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
430 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
431 function T = getBoundaryTractionOperator(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
432 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
433
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
434 T = obj.(['T', '_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
435 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
436
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
437 % 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
438 % 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
439 %
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
440 % boundary -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
441 function H = getBoundaryQuadrature(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
442 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 H = obj.getBoundaryQuadratureForScalarField(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
445 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
446 H = kron(H, I_dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
447 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
448
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
449 % 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
450 % 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
451 %
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
452 % boundary -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
453 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
454 assertIsMember(boundary, {'w', 'e', 's', 'n'})
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 H_b = obj.(['H_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
457 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
458
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
459 function N = size(obj)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
460 N = obj.dim*prod(obj.m);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well!
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
461 end
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