sbplib: +scheme/Elastic2dVariableAnisotropicUpwind.m annotate

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

Clean up hollow D2 generation by using precomputed coefficients stored in mat files

author	Martin Almquist <malmquist@stanford.edu>
date	Mon, 13 Jul 2020 14:21:41 -0700
parents	b5025bd67be1
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 Elastic2dVariableAnisotropicUpwind < scheme.Scheme
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	2
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	3 % Discretizes the elastic wave equation:
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	5 % opSet should be cell array of opSets, one per dimension. This
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	6 % is useful if we have periodic BC in one direction.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	7 % Assumes fully compatible operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	8
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	9 properties
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	10 m % Number of points in each direction, possibly a vector
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	11 h % Grid spacing
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	12
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	13 grid
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	14 dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	15
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	16 order % Order of accuracy for the approximation
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	17
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	18 % Diagonal matrices for variable coefficients
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	19 RHO, RHOi, RHOi_kron % Density
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	20 C % Elastic stiffness tensor
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	21
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	22 D % Total operator
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	23 Dp, Dm % First derivatives
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	24
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	25 % Boundary operators in cell format, used for BC
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	26 T_w, T_e, T_s, T_n
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	27
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	28 % Traction operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	29 tau_w, tau_e, tau_s, tau_n % Return vector field
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	30 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	31 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	32
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	33 % Inner products
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	34 H, Hi, Hi_kron, H_1D
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	35
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	36 % Boundary inner products (for scalar field)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	37 H_w, H_e, H_s, H_n
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	38
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	39 % Boundary restriction operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	40 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	41 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	42 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	43 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	44
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	45 % E{i}^T picks out component i
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	46 E
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	47
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	48 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	49 h11 % First entry in norm matrix
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	50
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	51 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	52
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	53 methods
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	54
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	55 % The coefficients can either be function handles or grid functions
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	56 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	57 function obj = Elastic2dVariableAnisotropicUpwind(g, order, rho, C, opSet, optFlag)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	58 default_arg('rho', @(x,y) 0*x+1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	59 default_arg('opSet',{@sbp.D1Upwind, @sbp.D1Upwind});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	60 default_arg('optFlag', false);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	61 dim = 2;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	62
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	63 C_default = cell(dim,dim,dim,dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	64 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	65 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	66 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	67 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	68 C_default{i,j,k,l} = @(x,y) 0*x + 1;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	69 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	70 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	71 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	72 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	73 default_arg('C', C_default);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	74 assert(isa(g, 'grid.Cartesian'))
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	75
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	76 if isa(rho, 'function_handle')
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	77 rho = grid.evalOn(g, rho);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	78 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	79
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	80 C_mat = cell(dim,dim,dim,dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	81 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	82 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	83 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	84 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	85 if isa(C{i,j,k,l}, 'function_handle')
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	86 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	87 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	88 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	89 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	90 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	91 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	92 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	93 obj.C = C_mat;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	94
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	95 m = g.size();
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	96 m_tot = g.N();
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	97 lim = g.lim;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	98 if isempty(lim)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	99 x = g.x;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	100 lim = cell(length(x),1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	101 for i = 1:length(x)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	102 lim{i} = {min(x{i}), max(x{i})};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	103 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	104 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	105
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	106 % 1D operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	107 ops = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	108 h = zeros(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	109 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	110 ops{i} = opSet{i}(m(i), lim{i}, order);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	111 h(i) = ops{i}.h;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	112 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	113
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	114 % Borrowing constants
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	115 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	116 obj.h11{i} = ops{i}.H(1,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	117 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	118
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	119 I = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	120 Dp = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	121 Dm = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	122 H = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	123 Hi = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	124 e_0 = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	125 e_m = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	126 d1_0 = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	127 d1_m = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	128
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	129 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	130 I{i} = speye(m(i));
1271 b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	131
b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	132 if isprop(ops{i}, 'Dp') && isprop(ops{i}, 'Dm')
b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	133 Dp{i} = ops{i}.Dp;
b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	134 Dm{i} = ops{i}.Dm;
b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	135 elseif isprop(ops{i}, 'D1')
b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	136 Dp{i} = ops{i}.D1;
b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	137 Dm{i} = ops{i}.D1;
b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	138 else
b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	139 error('opSet does not have Dp and Dm or D1');
b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	140 end
b5025bd67be1 Allow using D1 D1 in VariableAnisotropicUpwind Martin Almquist <malmquist@stanford.edu> parents: 1227 diff changeset	141
1227 02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	142 H{i} = ops{i}.H;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	143 Hi{i} = ops{i}.HI;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	144 e_0{i} = ops{i}.e_l;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	145 e_m{i} = ops{i}.e_r;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	146 d1_0{i} = (ops{i}.e_l' * Dm{i})';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	147 d1_m{i} = (ops{i}.e_r' * Dm{i})';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	148 end
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 %====== Assemble full operators ========
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	151 I_dim = speye(dim, dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	152 RHO = spdiag(rho);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	153 obj.RHO = RHO;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	154 obj.RHOi = inv(RHO);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	155 obj.RHOi_kron = kron(obj.RHOi, I_dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	156
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	157 obj.Dp = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	158 obj.Dm = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	159
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	160 % D1
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	161 obj.Dp{1} = kron(Dp{1},I{2});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	162 obj.Dp{2} = kron(I{1},Dp{2});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	163 obj.Dm{1} = kron(Dm{1},I{2});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	164 obj.Dm{2} = kron(I{1},Dm{2});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	165
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	166 % Boundary restriction operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	167 e_l = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	168 e_r = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	169 e_l{1} = kron(e_0{1}, I{2});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	170 e_l{2} = kron(I{1}, e_0{2});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	171 e_r{1} = kron(e_m{1}, I{2});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	172 e_r{2} = kron(I{1}, e_m{2});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	173
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	174 e_scalar_w = e_l{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	175 e_scalar_e = e_r{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	176 e_scalar_s = e_l{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	177 e_scalar_n = e_r{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	178
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	179 e_w = kron(e_scalar_w, I_dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	180 e_e = kron(e_scalar_e, I_dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	181 e_s = kron(e_scalar_s, I_dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	182 e_n = kron(e_scalar_n, I_dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	183
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	184 % 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	185 E = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	186 I = speye(m_tot,m_tot);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	187 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	188 e = sparse(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	189 e(i) = 1;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	190 E{i} = kron(I,e);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	191 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	192 obj.E = E;
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 e1_w = (e_scalar_w'*E{1}')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	195 e1_e = (e_scalar_e'*E{1}')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	196 e1_s = (e_scalar_s'*E{1}')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	197 e1_n = (e_scalar_n'*E{1}')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	198
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	199 e2_w = (e_scalar_w'*E{2}')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	200 e2_e = (e_scalar_e'*E{2}')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	201 e2_s = (e_scalar_s'*E{2}')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	202 e2_n = (e_scalar_n'*E{2}')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	203
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	204 % Quadratures
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	205 obj.H = kron(H{1},H{2});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	206 obj.Hi = inv(obj.H);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	207 obj.H_w = H{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	208 obj.H_e = H{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	209 obj.H_s = H{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	210 obj.H_n = H{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	211 obj.H_1D = {H{1}, H{2}};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	212
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	213 % Differentiation matrix D (without SAT)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	214 Dp = obj.Dp;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	215 Dm = obj.Dm;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	216 D = sparse(dimm_tot,dimm_tot);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	217 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	218 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	219 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	220 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	221 D = D + E{j}Dp{i}C_mat{i,j,k,l}Dm{k}E{l}';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	222 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	223 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	224 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	225 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	226 D = obj.RHOi_kron*D;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	227 obj.D = D;
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
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	230 % Numerical traction operators for BC.
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 % Formula at boundary j: % tau^{j}_i = sum_l T^{j}_{il} u_l
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	233 %
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	234 T_l = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	235 T_r = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	236 tau_l = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	237 tau_r = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	238
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	239 D1 = obj.Dm;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	240
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	241 % Boundary j
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	242 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	243 T_l{j} = cell(dim,dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	244 T_r{j} = cell(dim,dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	245 tau_l{j} = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	246 tau_r{j} = cell(dim,1);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	247
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	248 [~, n_l] = size(e_l{j});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	249 [~, n_r] = size(e_r{j});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	250
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	251 % Traction component i
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	252 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	253 tau_l{j}{i} = sparse(dim*m_tot, n_l);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	254 tau_r{j}{i} = sparse(dim*m_tot, n_r);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	255
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	256 % Displacement component l
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	257 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	258 T_l{j}{i,l} = sparse(m_tot, n_l);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	259 T_r{j}{i,l} = sparse(m_tot, n_r);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	260
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	261 % Derivative direction k
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	262 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	263 T_l{j}{i,l} = T_l{j}{i,l} ...
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	264 - (e_l{j}'C_mat{j,i,k,l}D1{k})';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	265 T_r{j}{i,l} = T_r{j}{i,l} ...
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	266 + (e_r{j}'C_mat{j,i,k,l}D1{k})';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	267 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	268 tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	269 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	270 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	271 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	272 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	273
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	274 % Traction tensors, T_ij
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	275 obj.T_w = T_l{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	276 obj.T_e = T_r{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	277 obj.T_s = T_l{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	278 obj.T_n = T_r{2};
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 % Restriction operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	281 obj.e_w = e_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	282 obj.e_e = e_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	283 obj.e_s = e_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	284 obj.e_n = e_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	285
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	286 obj.e1_w = e1_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	287 obj.e1_e = e1_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	288 obj.e1_s = e1_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	289 obj.e1_n = e1_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	290
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	291 obj.e2_w = e2_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	292 obj.e2_e = e2_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	293 obj.e2_s = e2_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	294 obj.e2_n = e2_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	295
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	296 obj.e_scalar_w = e_scalar_w;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	297 obj.e_scalar_e = e_scalar_e;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	298 obj.e_scalar_s = e_scalar_s;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	299 obj.e_scalar_n = e_scalar_n;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	300
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	301 % First component of traction
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	302 obj.tau1_w = tau_l{1}{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	303 obj.tau1_e = tau_r{1}{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	304 obj.tau1_s = tau_l{2}{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	305 obj.tau1_n = tau_r{2}{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	306
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	307 % Second component of traction
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	308 obj.tau2_w = tau_l{1}{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	309 obj.tau2_e = tau_r{1}{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	310 obj.tau2_s = tau_l{2}{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	311 obj.tau2_n = tau_r{2}{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	312
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	313 % Traction vectors
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	314 obj.tau_w = (e_w'e1_wobj.tau1_w')' + (e_w'e2_wobj.tau2_w')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	315 obj.tau_e = (e_e'e1_eobj.tau1_e')' + (e_e'e2_eobj.tau2_e')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	316 obj.tau_s = (e_s'e1_sobj.tau1_s')' + (e_s'e2_sobj.tau2_s')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	317 obj.tau_n = (e_n'e1_nobj.tau1_n')' + (e_n'e2_nobj.tau2_n')';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	318
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	319 % Kroneckered norms and coefficients
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	320 obj.Hi_kron = kron(obj.Hi, I_dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	321
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	322 % Misc.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	323 obj.m = m;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	324 obj.h = h;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	325 obj.order = order;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	326 obj.grid = g;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	327 obj.dim = dim;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	328
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	329 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	330
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	331
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	332 % 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	333 % 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	334 % 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	335 % 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	336 % 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	337 % {'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	338 % tangential/normal component.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	339 % 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	340 % 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	341 % 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	342
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	343 % For displacement bc:
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	344 % bc = {comp, 'd', dComps},
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	345 % where
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	346 % 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	347 % 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	348 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	349 default_arg('tuning', 1.0);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	350
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	351 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	352 comp = bc{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	353 type = bc{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	354 if ischar(comp)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	355 comp = obj.getComponent(comp, boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	356 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	357
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	358 e = obj.getBoundaryOperatorForScalarField('e', boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	359 tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	360 T = obj.getBoundaryTractionOperator(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	361 h11 = obj.getBorrowing(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	362 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	363 nu = obj.getNormal(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	364
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	365 E = obj.E;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	366 Hi = obj.Hi;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	367 RHOi = obj.RHOi;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	368 C = obj.C;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	369
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	370 dim = obj.dim;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	371 m_tot = obj.grid.N();
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	372
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	373 % Preallocate
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	374 [~, col] = size(tau);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	375 closure = sparse(dimm_tot, dimm_tot);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	376 penalty = sparse(dim*m_tot, col);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	377
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	378 j = comp;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	379 switch type
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 % Dirichlet boundary condition
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	382 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	383
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	384 if numel(bc) >= 3
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	385 dComps = bc{3};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	386 else
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	387 dComps = 1:dim;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	388 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	389
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	390 % Loops over components that Dirichlet penalties end up on
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	391 % Y: symmetrizing part of penalty
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	392 % Z: symmetric part of penalty
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	393 % X = Y + Z.
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 % Nonsymmetric part goes on all components to
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	396 % yield traction in discrete energy rate
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	397 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	398 Y = T{j,i}';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	399 X = e*Y;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	400 closure = closure + E{i}RHOiHiX'eH_gamma(e'*E{j}' );
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	401 penalty = penalty - E{i}RHOiHiX'e*H_gamma;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	402 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	403
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	404 % Symmetric part only required on components with displacement BC.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	405 % (Otherwise it's not symmetric.)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	406 for i = dComps
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	407 Z = sparse(m_tot, m_tot);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	408 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	409 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	410 Z = Z + nu(l)C{l,i,k,j}nu(k);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	411 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	412 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	413 Z = -tuningdim/h11Z;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	414 X = Z;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	415 closure = closure + E{i}RHOiHiX'eH_gamma(e'*E{j}' );
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	416 penalty = penalty - E{i}RHOiHiX'e*H_gamma;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	417 end
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 % Free boundary condition
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	420 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	421 closure = closure - E{j}RHOiHieH_gamma*tau';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	422 penalty = penalty + E{j}RHOiHieH_gamma;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	423
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	424 % Unknown boundary condition
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	425 otherwise
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	426 error('No such boundary condition: type = %s',type);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	427 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	428 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	429
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	430 % 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	431 % Fields:
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	432 % -- 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	433 % -- interpolation: type of interpolation, default 'none'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	434 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	435
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	436 defaultType.tuning = 1.0;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	437 defaultType.interpolation = 'none';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	438 default_struct('type', defaultType);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	439
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	440 switch type.interpolation
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	441 case {'none', ''}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	442 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	443 case {'op','OP'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	444 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	445 otherwise
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	446 error('Unknown type of interpolation: %s ', type.interpolation);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	447 end
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 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	451 tuning = type.tuning;
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 % u denotes the solution in the own domain
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	454 % v denotes the solution in the neighbour domain
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	455
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	456 u = obj;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	457 v = neighbour_scheme;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	458
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	459 % Operators, u side
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	460 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	461 tau_u = u.getBoundaryOperator('tau', boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	462 h11_u = u.getBorrowing(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	463 nu_u = u.getNormal(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	464
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	465 E_u = u.E;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	466 C_u = u.C;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	467 m_tot_u = u.grid.N();
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	468
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	469 % Operators, v side
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	470 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	471 tau_v = v.getBoundaryOperator('tau', neighbour_boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	472 h11_v = v.getBorrowing(neighbour_boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	473 nu_v = v.getNormal(neighbour_boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	474
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	475 E_v = v.E;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	476 C_v = v.C;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	477 m_tot_v = v.grid.N();
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	478
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	479 % Operators that are only required for own domain
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	480 Hi = u.Hi_kron;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	481 RHOi = u.RHOi_kron;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	482 e_kron = u.getBoundaryOperator('e', boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	483 T_u = u.getBoundaryTractionOperator(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	484
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	485 % Shared operators
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	486 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	487 H_gamma_kron = u.getBoundaryQuadrature(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	488 dim = u.dim;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	489
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	490 % Preallocate
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	491 [~, m_int] = size(H_gamma);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	492 closure = sparse(dimm_tot_u, dimm_tot_u);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	493 penalty = sparse(dimm_tot_u, dimm_tot_v);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	494
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	495 % ---- Continuity of displacement ------
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	496
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	497 % Y: symmetrizing part of penalty
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	498 % Z: symmetric part of penalty
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	499 % X = Y + Z.
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	500
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	501 % Loop over components to couple across interface
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	502 for j = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	503
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	504 % Loop over components that penalties end up on
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	505 for i = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	506 Y = 1/2*T_u{j,i}';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	507 Z_u = sparse(m_int, m_int);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	508 Z_v = sparse(m_int, m_int);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	509 for l = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	510 for k = 1:dim
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	511 Z_u = Z_u + e_u'nu_u(l)C_u{l,i,k,j}nu_u(k)e_u;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	512 Z_v = Z_v + e_v'nu_v(l)C_v{l,i,k,j}nu_v(k)e_v;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	513 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	514 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	515 Z = -tuningdim( 1/(4h11_u)Z_u + 1/(4h11_v)Z_v );
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	516 X = Y + Z*e_u';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	517 closure = closure + E_u{i}X'H_gammae_u'E_u{j}';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	518 penalty = penalty - E_u{i}X'H_gammae_v'E_v{j}';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	519 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	520 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	521
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	522 % ---- Continuity of traction ------
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	523 closure = closure - 1/2e_kronH_gamma_kron*tau_u';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	524 penalty = penalty - 1/2e_kronH_gamma_kron*tau_v';
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	525
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	526 % ---- Multiply by inverse of density x quadraure ----
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	527 closure = RHOiHiclosure;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	528 penalty = RHOiHipenalty;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	529
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	530 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	531
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	532 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	533 error('Non-conforming interfaces not implemented yet.');
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	534 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	535
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	536 % Returns the component number that is the tangential/normal component
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	537 % at the specified boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	538 function comp = getComponent(obj, comp_str, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	539 assertIsMember(comp_str, {'normal', 'tangential'});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	540 assertIsMember(boundary, {'w', 'e', 's', 'n'});
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	541
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	542 switch boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	543 case {'w', 'e'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	544 switch comp_str
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	545 case 'normal'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	546 comp = 1;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	547 case 'tangential'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	548 comp = 2;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	549 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	550 case {'s', 'n'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	551 switch comp_str
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	552 case 'normal'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	553 comp = 2;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	554 case 'tangential'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	555 comp = 1;
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	556 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	557 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	558 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	559
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	560 % 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	561 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	562 function h11 = getBorrowing(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	563 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	564
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	565 switch boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	566 case {'w','e'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	567 h11 = obj.h11{1};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	568 case {'s', 'n'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	569 h11 = obj.h11{2};
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	570 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	571 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	572
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	573 % 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	574 function nu = getNormal(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	575 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	576
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	577 switch boundary
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	578 case 'w'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	579 nu = [-1,0];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	580 case 'e'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	581 nu = [1,0];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	582 case 's'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	583 nu = [0,-1];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	584 case 'n'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	585 nu = [0,1];
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	586 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	587 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	588
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	589 % 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	590 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	591 function o = getBoundaryOperator(obj, op, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	592 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	593 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	594
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	595 switch op
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	596 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	597 o = obj.([op, '_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	598 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	599
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	600 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	601
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	602 % 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	603 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	604 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	605 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	606 assertIsMember(op, {'e'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	607
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	608 switch op
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	609
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	610 case 'e'
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	611 o = obj.(['e_scalar', '_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	612 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	613
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	614 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	615
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	616 % Returns 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	617 % 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	618 % op -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	619 function T = getBoundaryTractionOperator(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	620 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	621
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	622 T = obj.(['T', '_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	623 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	624
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	625 % 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	626 % 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	627 %
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	628 % boundary -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	629 function H = getBoundaryQuadrature(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	630 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	631
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	632 H = obj.getBoundaryQuadratureForScalarField(boundary);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	633 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	634 H = kron(H, I_dim);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	635 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	636
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	637 % 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	638 % 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	639 %
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	640 % boundary -- string
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	641 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	642 assertIsMember(boundary, {'w', 'e', 's', 'n'})
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	643
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	644 H_b = obj.(['H_', boundary]);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	645 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	646
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	647 function N = size(obj)
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	648 N = obj.dim*prod(obj.m);
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	649 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	650 end
02dfe3a56742 Add Upwind ElasticAnisotropic schemes. Seem to work really well! Martin Almquist <malmquist@stanford.edu> parents: diff changeset	651 end

Mercurial > repos > public > sbplib

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