sbplib: +scheme/Elastic2dCurvilinearAnisotropic.m annotate

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1295:cb053fabbedc feature/poroelastic

Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component.

author	Martin Almquist <malmquist@stanford.edu>
date	Thu, 02 Jul 2020 20:05:50 -0700
parents	fe712a1fca3f
children	a0d615bde7f8

rev	line source
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	1 classdef Elastic2dCurvilinearAnisotropic < scheme.Scheme
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	2
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	3 % Discretizes the elastic wave equation:
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	5 % in curvilinear coordinates.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	6 % opSet should be cell array of opSets, one per dimension. This
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	7 % is useful if we have periodic BC in one direction.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	8 % Assumes fully compatible operators.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	9
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	10 properties
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	11 m % Number of points in each direction, possibly a vector
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	12 h % Grid spacing
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	13
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	14 grid
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	15 dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	16
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	17 order % Order of accuracy for the approximation
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	18
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	19 % Diagonal matrices for variable coefficients
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	20 J, Ji
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	21 RHO % Density
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	22 C % Elastic stiffness tensor
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	23
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	24 D % Total operator
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	25
1211 3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1209 diff changeset	26 Dx, Dy % Physical derivatives
1214 e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	27 sigma % Cell matrix of physical stress operators
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	28 n_w, n_e, n_s, n_n % Physical normals
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	29 tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents
1211 3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1209 diff changeset	30
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	31 % Boundary operators in cell format, used for BC
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	32 T_w, T_e, T_s, T_n
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	33
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	34 % Traction operators
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	35 tau_w, tau_e, tau_s, tau_n % Return vector field
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	36 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	37 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	38 tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	39 tau_t_w, tau_t_e, tau_t_s, tau_t_n % Return scalar field
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	40
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	41 % Inner products
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	42 H
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	43
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	44 % Boundary inner products (for scalar field)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	45 H_w, H_e, H_s, H_n
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	46
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	47 % Surface Jacobian vectors
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	48 s_w, s_e, s_s, s_n
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	49
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	50 % Boundary restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	51 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	52 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	53 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	54 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	55 en_w, en_e, en_s, en_n % Act on vector field, return normal component
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	56 et_w, et_e, et_s, et_n % Act on vector field, return tangential component
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	57
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	58 % E{i}^T picks out component i
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	59 E
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	60
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	61 % Elastic2dVariableAnisotropic object for reference domain
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	62 refObj
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	63 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	64
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	65 methods
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	66
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	67 % The coefficients can either be function handles or grid functions
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	68 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	69 function obj = Elastic2dCurvilinearAnisotropic(g, order, rho, C, opSet, optFlag)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	70 default_arg('rho', @(x,y) 0*x+1);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	71 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	72 default_arg('optFlag', false);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	73 dim = 2;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	74
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	75 C_default = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	76 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	77 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	78 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	79 for l = 1:dim
1272 15865fbda16e Set default stiffness tensors to 0 in CurvilinearAnisotropic schemes Martin Almquist <malmquist@stanford.edu> parents: 1214 diff changeset	80 C_default{i,j,k,l} = @(x,y) 0*x ;
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	81 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	82 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	83 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	84 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	85 default_arg('C', C_default);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	86
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	87 assert(isa(g, 'grid.Curvilinear'));
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	88
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	89 if isa(rho, 'function_handle')
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	90 rho = grid.evalOn(g, rho);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	91 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	92
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	93 C_mat = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	94 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	95 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	96 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	97 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	98 if isa(C{i,j,k,l}, 'function_handle')
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	99 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	100 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	101 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	102 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	103 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	104 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	105 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	106 obj.C = C_mat;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	107
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	108 m = g.size();
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	109 m_tot = g.N();
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	110
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	111 % 1D operators
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	112 m_u = m(1);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	113 m_v = m(2);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	114 ops_u = opSet{1}(m_u, {0, 1}, order);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	115 ops_v = opSet{2}(m_v, {0, 1}, order);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	116
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	117 h_u = ops_u.h;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	118 h_v = ops_v.h;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	119
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	120 I_u = speye(m_u);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	121 I_v = speye(m_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	122
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	123 D1_u = ops_u.D1;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	124 H_u = ops_u.H;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	125 Hi_u = ops_u.HI;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	126 e_l_u = ops_u.e_l;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	127 e_r_u = ops_u.e_r;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	128 d1_l_u = ops_u.d1_l;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	129 d1_r_u = ops_u.d1_r;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	130
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	131 D1_v = ops_v.D1;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	132 H_v = ops_v.H;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	133 Hi_v = ops_v.HI;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	134 e_l_v = ops_v.e_l;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	135 e_r_v = ops_v.e_r;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	136 d1_l_v = ops_v.d1_l;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	137 d1_r_v = ops_v.d1_r;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	138
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	139
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	140 % Logical operators
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	141 Du = kr(D1_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	142 Dv = kr(I_u,D1_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	143
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	144 e_w = kr(e_l_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	145 e_e = kr(e_r_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	146 e_s = kr(I_u,e_l_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	147 e_n = kr(I_u,e_r_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	148
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	149 % Metric coefficients
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	150 coords = g.points();
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	151 x = coords(:,1);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	152 y = coords(:,2);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	153
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	154 x_u = Du*x;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	155 x_v = Dv*x;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	156 y_u = Du*y;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	157 y_v = Dv*y;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	158
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	159 J = x_u.y_v - x_v.y_u;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	160
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	161 K = cell(dim, dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	162 K{1,1} = y_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	163 K{1,2} = -y_u./J;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	164 K{2,1} = -x_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	165 K{2,2} = x_u./J;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	166
1211 3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1209 diff changeset	167 % Physical derivatives
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1209 diff changeset	168 obj.Dx = spdiag( y_v./J)Du + spdiag(-y_u./J)Dv;
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1209 diff changeset	169 obj.Dy = spdiag(-x_v./J)Du + spdiag( x_u./J)Dv;
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1209 diff changeset	170
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	171 % Wrap around Aniosotropic Cartesian
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	172 rho_tilde = J.*rho;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	173
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	174 PHI = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	175 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	176 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	177 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	178 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	179 PHI{i,j,k,l} = 0*C{i,j,k,l};
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	180 for m = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	181 for n = 1:dim
1209 67eee83fd9c9 Swap indices in anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1205 diff changeset	182 PHI{i,j,k,l} = PHI{i,j,k,l} + J.K{m,i}.C{m,j,n,l}.*K{n,k};
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	183 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	184 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	185 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	186 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	187 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	188 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	189
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	190 gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1});
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	191 refObj = scheme.Elastic2dVariableAnisotropic(gRef, order, rho_tilde, PHI, opSet);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	192
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	193 %---- Set object properties ------
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	194 obj.RHO = spdiag(rho);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	195
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	196 % Volume quadrature
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	197 obj.J = spdiag(J);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	198 obj.Ji = spdiag(1./J);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	199 obj.H = obj.J*kr(H_u,H_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	200
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	201 % Boundary quadratures
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	202 s_w = sqrt((e_w'x_v).^2 + (e_w'y_v).^2);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	203 s_e = sqrt((e_e'x_v).^2 + (e_e'y_v).^2);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	204 s_s = sqrt((e_s'x_u).^2 + (e_s'y_u).^2);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	205 s_n = sqrt((e_n'x_u).^2 + (e_n'y_u).^2);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	206 obj.s_w = s_w;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	207 obj.s_e = s_e;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	208 obj.s_s = s_s;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	209 obj.s_n = s_n;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	210
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	211 obj.H_w = H_v*spdiag(s_w);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	212 obj.H_e = H_v*spdiag(s_e);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	213 obj.H_s = H_u*spdiag(s_s);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	214 obj.H_n = H_u*spdiag(s_n);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	215
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	216 % Restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	217 obj.e_w = refObj.e_w;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	218 obj.e_e = refObj.e_e;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	219 obj.e_s = refObj.e_s;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	220 obj.e_n = refObj.e_n;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	221
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	222 % Adapt things from reference object
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	223 obj.D = refObj.D;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	224 obj.E = refObj.E;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	225
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	226 obj.e1_w = refObj.e1_w;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	227 obj.e1_e = refObj.e1_e;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	228 obj.e1_s = refObj.e1_s;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	229 obj.e1_n = refObj.e1_n;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	230
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	231 obj.e2_w = refObj.e2_w;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	232 obj.e2_e = refObj.e2_e;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	233 obj.e2_s = refObj.e2_s;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	234 obj.e2_n = refObj.e2_n;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	235
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	236 obj.e_scalar_w = refObj.e_scalar_w;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	237 obj.e_scalar_e = refObj.e_scalar_e;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	238 obj.e_scalar_s = refObj.e_scalar_s;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	239 obj.e_scalar_n = refObj.e_scalar_n;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	240
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	241 e1_w = obj.e1_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	242 e1_e = obj.e1_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	243 e1_s = obj.e1_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	244 e1_n = obj.e1_n;
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	245
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	246 e2_w = obj.e2_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	247 e2_e = obj.e2_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	248 e2_s = obj.e2_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	249 e2_n = obj.e2_n;
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	250
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	251 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	252 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	253 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	254 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	255
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	256 obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	257 obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	258 obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	259 obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	260
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	261 obj.tau_w = (refObj.e_w'obj.e1_wobj.tau1_w')' + (refObj.e_w'obj.e2_wobj.tau2_w')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	262 obj.tau_e = (refObj.e_e'obj.e1_eobj.tau1_e')' + (refObj.e_e'obj.e2_eobj.tau2_e')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	263 obj.tau_s = (refObj.e_s'obj.e1_sobj.tau1_s')' + (refObj.e_s'obj.e2_sobj.tau2_s')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	264 obj.tau_n = (refObj.e_n'obj.e1_nobj.tau1_n')' + (refObj.e_n'obj.e2_nobj.tau2_n')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	265
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	266 % Physical normals
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	267 e_w = obj.e_scalar_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	268 e_e = obj.e_scalar_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	269 e_s = obj.e_scalar_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	270 e_n = obj.e_scalar_n;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	271
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	272 e_w_vec = obj.e_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	273 e_e_vec = obj.e_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	274 e_s_vec = obj.e_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	275 e_n_vec = obj.e_n;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	276
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	277 nu_w = [-1,0];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	278 nu_e = [1,0];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	279 nu_s = [0,-1];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	280 nu_n = [0,1];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	281
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	282 obj.n_w = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	283 obj.n_e = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	284 obj.n_s = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	285 obj.n_n = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	286
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	287 % Compute normal and rotate (exactly!) 90 degrees counter-clockwise to get tangent
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	288 n_w_1 = (1./s_w).e_w'(J.(K{1,1}nu_w(1) + K{1,2}*nu_w(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	289 n_w_2 = (1./s_w).e_w'(J.(K{2,1}nu_w(1) + K{2,2}*nu_w(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	290 obj.n_w{1} = spdiag(n_w_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	291 obj.n_w{2} = spdiag(n_w_2);
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	292 obj.tangent_w = {-obj.n_w{2}, obj.n_w{1}};
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	293
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	294 n_e_1 = (1./s_e).e_e'(J.(K{1,1}nu_e(1) + K{1,2}*nu_e(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	295 n_e_2 = (1./s_e).e_e'(J.(K{2,1}nu_e(1) + K{2,2}*nu_e(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	296 obj.n_e{1} = spdiag(n_e_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	297 obj.n_e{2} = spdiag(n_e_2);
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	298 obj.tangent_e = {-obj.n_e{2}, obj.n_e{1}};
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	299
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	300 n_s_1 = (1./s_s).e_s'(J.(K{1,1}nu_s(1) + K{1,2}*nu_s(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	301 n_s_2 = (1./s_s).e_s'(J.(K{2,1}nu_s(1) + K{2,2}*nu_s(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	302 obj.n_s{1} = spdiag(n_s_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	303 obj.n_s{2} = spdiag(n_s_2);
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	304 obj.tangent_s = {-obj.n_s{2}, obj.n_s{1}};
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	305
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	306 n_n_1 = (1./s_n).e_n'(J.(K{1,1}nu_n(1) + K{1,2}*nu_n(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	307 n_n_2 = (1./s_n).e_n'(J.(K{2,1}nu_n(1) + K{2,2}*nu_n(2)));
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	308 obj.n_n{1} = spdiag(n_n_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	309 obj.n_n{2} = spdiag(n_n_2);
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	310 obj.tangent_n = {-obj.n_n{2}, obj.n_n{1}};
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	311
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	312 % Operators that extract the normal component
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	313 obj.en_w = (obj.n_w{1}obj.e1_w' + obj.n_w{2}obj.e2_w')';
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	314 obj.en_e = (obj.n_e{1}obj.e1_e' + obj.n_e{2}obj.e2_e')';
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	315 obj.en_s = (obj.n_s{1}obj.e1_s' + obj.n_s{2}obj.e2_s')';
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	316 obj.en_n = (obj.n_n{1}obj.e1_n' + obj.n_n{2}obj.e2_n')';
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	317
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	318 % Operators that extract the tangential component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	319 obj.et_w = (obj.tangent_w{1}obj.e1_w' + obj.tangent_w{2}obj.e2_w')';
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	320 obj.et_e = (obj.tangent_e{1}obj.e1_e' + obj.tangent_e{2}obj.e2_e')';
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	321 obj.et_s = (obj.tangent_s{1}obj.e1_s' + obj.tangent_s{2}obj.e2_s')';
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	322 obj.et_n = (obj.tangent_n{1}obj.e1_n' + obj.tangent_n{2}obj.e2_n')';
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	323
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	324 obj.tau_n_w = (obj.n_w{1}obj.tau1_w' + obj.n_w{2}obj.tau2_w')';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	325 obj.tau_n_e = (obj.n_e{1}obj.tau1_e' + obj.n_e{2}obj.tau2_e')';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	326 obj.tau_n_s = (obj.n_s{1}obj.tau1_s' + obj.n_s{2}obj.tau2_s')';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	327 obj.tau_n_n = (obj.n_n{1}obj.tau1_n' + obj.n_n{2}obj.tau2_n')';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	328
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	329 obj.tau_t_w = (obj.tangent_w{1}obj.tau1_w' + obj.tangent_w{2}obj.tau2_w')';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	330 obj.tau_t_e = (obj.tangent_e{1}obj.tau1_e' + obj.tangent_e{2}obj.tau2_e')';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	331 obj.tau_t_s = (obj.tangent_s{1}obj.tau1_s' + obj.tangent_s{2}obj.tau2_s')';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	332 obj.tau_t_n = (obj.tangent_n{1}obj.tau1_n' + obj.tangent_n{2}obj.tau2_n')';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	333
1214 e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	334 % Stress operators
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	335 sigma = cell(dim, dim);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	336 D1 = {obj.Dx, obj.Dy};
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	337 E = obj.E;
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	338 N = length(obj.RHO);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	339 for i = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	340 for j = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	341 sigma{i,j} = sparse(N,2*N);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	342 for k = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	343 for l = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	344 sigma{i,j} = sigma{i,j} + obj.C{i,j,k,l}D1{k}E{l}';
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	345 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	346 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	347 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	348 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	349 obj.sigma = sigma;
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	350
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	351 % Misc.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	352 obj.refObj = refObj;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	353 obj.m = refObj.m;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	354 obj.h = refObj.h;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	355 obj.order = order;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	356 obj.grid = g;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	357 obj.dim = dim;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	358
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	359 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	360
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	361
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	362 % Closure functions return the operators applied to the own domain to close the boundary
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	363 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	364 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	365 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	366 % on the first component. Can also be e.g.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	367 % {'normal', 'd'} or {'tangential', 't'} for conditions on
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	368 % tangential/normal component.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	369 % data is a function returning the data that should be applied at the boundary.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	370 % neighbour_scheme is an instance of Scheme that should be interfaced to.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	371 % neighbour_boundary is a string specifying which boundary to interface to.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	372
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	373 % For displacement bc:
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	374 % bc = {comp, 'd', dComps},
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	375 % where
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	376 % dComps = vector of components with displacement BC. Default: 1:dim.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	377 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	378 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	379 default_arg('tuning', 1.0);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	380 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	381
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	382 component = bc{1};
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	383 type = bc{2};
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	384
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	385 switch component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	386
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	387 % If conditions on Cartesian components
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	388 case {1,2}
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	389 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	390
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	391 % If conditions on normal or tangential components
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	392 case {'n', 't'}
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	393
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	394 switch component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	395 case 'n'
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	396 en = obj.getBoundaryOperator('en', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	397 case 't'
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	398 en = obj.getBoundaryOperator('et', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	399 end
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	400 e1 = obj.getBoundaryOperator('e1', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	401
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	402 bc1 = {1, type};
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	403 [c1, p1] = obj.refObj.boundary_condition(boundary, bc1, tuning);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	404 bc2 = {2, type};
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	405 c2 = obj.refObj.boundary_condition(boundary, bc2, tuning);
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	406
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	407 switch type
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	408 case {'F','f','Free','free','traction','Traction','t','T'}
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	409 closure = enen'(c1+c2);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	410 penalty = ene1'p1;
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	411 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	412 [closure, penalty] = obj.displacementBCNormalTangential(boundary, bc, tuning);
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	413 end
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	414
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	415 end
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	416
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	417 switch type
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	418 case {'F','f','Free','free','traction','Traction','t','T'}
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	419
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	420 s = obj.(['s_' boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	421 s = spdiag(s);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	422 penalty = penalty*s;
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	423
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	424 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	425 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	426
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	427 function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning)
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	428 disp('WARNING: displacementBCNormalTangential is only guaranteed to work for displaacement BC on one component and traction on the other');
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	429 u = obj;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	430
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	431 component = bc{1};
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	432 type = bc{2};
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	433
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	434 switch component
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	435 case 'n'
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	436 en = u.getBoundaryOperator('en', boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	437 tau_n = u.getBoundaryOperator('tau_n', boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	438 N = u.getNormal(boundary);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	439 case 't'
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	440 en = u.getBoundaryOperator('et', boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	441 tau_n = u.getBoundaryOperator('tau_t', boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	442 N = u.getTangent(boundary);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	443 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	444
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	445 % Operators
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	446 e = u.getBoundaryOperatorForScalarField('e', boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	447 h11 = u.getBorrowing(boundary);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	448 n = u.getNormal(boundary);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	449
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	450 C = u.C;
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	451 Ji = u.Ji;
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	452 s = spdiag(u.(['s_' boundary]));
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	453 m_tot = u.grid.N();
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	454
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	455 Hi = u.E{1}inv(u.H)u.E{1}' + u.E{2}inv(u.H)u.E{2}';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	456 RHOi = u.E{1}inv(u.RHO)u.E{1}' + u.E{2}inv(u.RHO)u.E{2}';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	457
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	458 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	459 dim = u.dim;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	460
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	461 % Preallocate
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	462 [~, m_int] = size(H_gamma);
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	463 closure = sparse(dimm_tot, dimm_tot);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	464 penalty = sparse(dim*m_tot, m_int);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	465
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	466 % Term 1: The symmetric term
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	467 Z = sparse(m_int, m_int);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	468 for i = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	469 for j = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	470 for l = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	471 for k = 1:dim
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	472 Z = Z + n{i}N{j}e'JiC{i,j,k,l}en{k}*N{l};
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	473 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	474 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	475 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	476 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	477
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	478 Z = -tuningdim1/h11sZ;
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	479 closure = closure + enH_gammaZ*en';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	480 penalty = penalty - enH_gammaZ;
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	481
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	482 % Term 2: The symmetrizing term
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	483 closure = closure + tau_nH_gammaen';
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	484 penalty = penalty - tau_n*H_gamma;
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	485
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	486 % Multiply all normal component terms by inverse of density x quadraure
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	487 closure = RHOiHiclosure;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	488 penalty = RHOiHipenalty;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	489 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	490
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	491 % type Struct that specifies the interface coupling.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	492 % Fields:
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	493 % -- tuning: penalty strength, defaults to 1.0
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	494 % -- interpolation: type of interpolation, default 'none'
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	495 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	496
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	497 defaultType.tuning = 1.0;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	498 defaultType.interpolation = 'none';
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	499 defaultType.type = 'standard';
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	500 default_struct('type', defaultType);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	501
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	502 switch type.type
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	503 case 'standard'
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	504 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	505 case 'frictionalFault'
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	506 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	507 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	508
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	509 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	510
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	511 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	512 tuning = type.tuning;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	513
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	514 % u denotes the solution in the own domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	515 % v denotes the solution in the neighbour domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	516
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	517 u = obj;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	518 v = neighbour_scheme;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	519
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	520 % Operators, u side
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	521 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	522 en_u = u.getBoundaryOperator('en', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	523 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	524 h11_u = u.getBorrowing(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	525 n_u = u.getNormal(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	526
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	527 C_u = u.C;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	528 Ji_u = u.Ji;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	529 s_u = spdiag(u.(['s_' boundary]));
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	530 m_tot_u = u.grid.N();
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	531
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	532 % Operators, v side
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	533 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	534 en_v = v.getBoundaryOperator('en', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	535 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	536 h11_v = v.getBorrowing(neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	537 n_v = v.getNormal(neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	538
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	539 C_v = v.C;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	540 Ji_v = v.Ji;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	541 s_v = spdiag(v.(['s_' neighbour_boundary]));
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	542 m_tot_v = v.grid.N();
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	543
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	544 % Operators that are only required for own domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	545 Hi = u.E{1}inv(u.H)u.E{1}' + u.E{2}inv(u.H)u.E{2}';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	546 RHOi = u.E{1}inv(u.RHO)u.E{1}' + u.E{2}inv(u.RHO)u.E{2}';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	547
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	548 % Shared operators
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	549 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	550 dim = u.dim;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	551
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	552 % Preallocate
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	553 [~, m_int] = size(H_gamma);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	554 closure = sparse(dimm_tot_u, dimm_tot_u);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	555 penalty = sparse(dimm_tot_u, dimm_tot_v);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	556
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	557 % Continuity of normal displacement, term 1: The symmetric term
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	558 Z_u = sparse(m_int, m_int);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	559 Z_v = sparse(m_int, m_int);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	560 for i = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	561 for j = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	562 for l = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	563 for k = 1:dim
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	564 Z_u = Z_u + n_u{i}n_u{j}e_u'Ji_uC_u{i,j,k,l}e_un_u{k}*n_u{l};
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	565 Z_v = Z_v + n_v{i}n_v{j}e_v'Ji_vC_v{i,j,k,l}e_vn_v{k}*n_v{l};
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	566 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	567 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	568 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	569 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	570
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	571 Z = -tuningdim( 1/(4h11_u)s_uZ_u + 1/(4h11_v)s_vZ_v );
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	572 closure = closure + en_uH_gammaZ*en_u';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	573 penalty = penalty + en_uH_gammaZ*en_v';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	574
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	575 % Continuity of normal displacement, term 2: The symmetrizing term
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	576 closure = closure + 1/2tau_n_uH_gamma*en_u';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	577 penalty = penalty + 1/2tau_n_uH_gamma*en_v';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	578
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	579 % Continuity of normal traction
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	580 closure = closure - 1/2en_uH_gamma*tau_n_u';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	581 penalty = penalty + 1/2en_uH_gamma*tau_n_v';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	582
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	583 % Multiply all normal component terms by inverse of density x quadraure
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	584 closure = RHOiHiclosure;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	585 penalty = RHOiHipenalty;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	586
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	587 % ---- Tangential tractions are imposed just like traction BC ------
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	588 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	589
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	590 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	591
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	592
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	593 % Returns h11 for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	594 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	595 function h11 = getBorrowing(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	596 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	597
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	598 switch boundary
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	599 case {'w','e'}
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	600 h11 = obj.refObj.h11{1};
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	601 case {'s', 'n'}
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	602 h11 = obj.refObj.h11{2};
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	603 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	604 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	605
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	606 % Returns the outward unit normal vector for the boundary specified by the string boundary.
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	607 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	608 function n = getNormal(obj, boundary)
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	609 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	610
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	611 n = obj.(['n_' boundary]);
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	612 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	613
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	614 % Returns the unit tangent vector for the boundary specified by the string boundary.
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	615 % t is a cell of diagonal matrices for each normal component, t{1} = t_1, t{2} = t_2.
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	616 function t = getTangent(obj, boundary)
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	617 assertIsMember(boundary, {'w', 'e', 's', 'n'})
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	618
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	619 t = obj.(['tangent_' boundary]);
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	620 end
cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	621
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	622 % Returns the boundary operator op for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	623 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	624 function o = getBoundaryOperator(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	625 assertIsMember(boundary, {'w', 'e', 's', 'n'})
1295 cb053fabbedc Add displacement BC on normal tangential components, but it is only provably stable if traction bc are used for the other component. Martin Almquist <malmquist@stanford.edu> parents: 1294 diff changeset	626 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	627
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	628 o = obj.([op, '_', boundary]);
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	629
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	630 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	631
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	632 % Returns the boundary operator op for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	633 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	634 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	635 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	636 assertIsMember(op, {'e'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	637
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	638 switch op
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	639
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	640 case 'e'
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	641 o = obj.(['e_scalar', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	642 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	643
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	644 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	645
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	646 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	647 % Formula: tau_i = T_ij u_j
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	648 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	649 function T = getBoundaryTractionOperator(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	650 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	651
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	652 T = obj.(['T', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	653 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	654
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	655 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	656 % corresponding to the number of boundary unknowns
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	657 %
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	658 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	659 function H = getBoundaryQuadrature(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	660 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	661
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	662 H = obj.getBoundaryQuadratureForScalarField(boundary);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	663 I_dim = speye(obj.dim, obj.dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	664 H = kron(H, I_dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	665 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	666
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	667 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	668 % corresponding to the number of boundary grid points
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	669 %
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	670 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	671 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	672 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	673
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	674 H_b = obj.(['H_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	675 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	676
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	677 function N = size(obj)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	678 N = obj.dim*prod(obj.m);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	679 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	680 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	681 end

Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1295:cb053fabbedc feature/poroelastic