sbplib: +scheme/Elastic2dCurvilinearAnisotropicUpwind.m annotate

annotate +scheme/Elastic2dCurvilinearAnisotropicUpwind.m @ 1260:7fbee3de0ab2 feature/poroelastic

Remove dummy line in Eastic2dVariableAnisotropic.

author	Martin Almquist <malmquist@stanford.edu>
date	Wed, 29 Apr 2020 21:01:09 -0700
parents	02dfe3a56742
children	9a7cac4202d0

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

Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropicUpwind.m @ 1260:7fbee3de0ab2 feature/poroelastic