sbplib: +scheme/Elastic2dCurvilinearAnisotropicUpwind.m annotate

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

Fix comment

author	Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date	Fri, 22 Jul 2022 16:31:18 +0200
parents	8ff3a95ad7cc
children

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_wobj.tau1_w')' + (refObj.e_w'obj.e2_wobj.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_eobj.tau1_e')' + (refObj.e_e'obj.e2_eobj.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_sobj.tau1_s')' + (refObj.e_s'obj.e2_sobj.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_nobj.tau1_n')' + (refObj.e_n'obj.e2_nobj.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

Mercurial > repos > public > sbplib

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