sbplib: +scheme/Elastic2dCurvilinearAnisotropic.m annotate

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1319:8b1110385ee2 feature/poroelastic

Bugfix interface forcing in ElasticCurveAnisotropic

author	Martin Almquist <malmquist@stanford.edu>
date	Tue, 28 Jul 2020 21:58:22 -0700
parents	34997aced843
children	b40359c9faed

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
1306 633757e582e5 Add transformation gradient as property in CurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	26 K % Transformation gradient
1211 3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1209 diff changeset	27 Dx, Dy % Physical derivatives
1214 e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	28 sigma % Cell matrix of physical stress operators
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	29 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	30 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	31
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	32 % Boundary operators in cell format, used for BC
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	33 T_w, T_e, T_s, T_n
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	34
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	35 % Traction operators
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	36 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	37 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	38 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	39 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	40 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	41
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	42 % Inner products
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	43 H
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	44
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	45 % Boundary inner products (for scalar field)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	46 H_w, H_e, H_s, H_n
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	47
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	48 % Surface Jacobian vectors
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	49 s_w, s_e, s_s, s_n
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	50
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	51 % Boundary restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	52 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	53 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	54 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	55 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	56 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	57 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	58
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	59 % E{i}^T picks out component i
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	60 E
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	61
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	62 % Elastic2dVariableAnisotropic object for reference domain
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	63 refObj
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	64 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	65
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	66 methods
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	67
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	68 % The coefficients can either be function handles or grid functions
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	69 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1295 diff changeset	70 function obj = Elastic2dCurvilinearAnisotropic(g, order, rho, C, opSet, optFlag, hollow)
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1295 diff changeset	71 default_arg('hollow', false);
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	72 default_arg('rho', @(x,y) 0*x+1);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	73 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	74 default_arg('optFlag', false);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	75 dim = 2;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	76
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	77 C_default = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	78 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	79 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	80 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	81 for l = 1:dim
1272 15865fbda16e Set default stiffness tensors to 0 in CurvilinearAnisotropic schemes Martin Almquist <malmquist@stanford.edu> parents: 1214 diff changeset	82 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	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 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	86 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	87 default_arg('C', C_default);
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 assert(isa(g, 'grid.Curvilinear'));
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	90
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	91 if isa(rho, 'function_handle')
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	92 rho = grid.evalOn(g, rho);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	93 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	94
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	95 C_mat = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	96 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	97 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	98 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	99 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	100 if isa(C{i,j,k,l}, 'function_handle')
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	101 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	102 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	103 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	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 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	107 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	108 obj.C = C_mat;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	109
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	110 m = g.size();
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	111 m_tot = g.N();
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	112
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	113 % 1D operators
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	114 m_u = m(1);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	115 m_v = m(2);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	116 ops_u = opSet{1}(m_u, {0, 1}, order);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	117 ops_v = opSet{2}(m_v, {0, 1}, order);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	118
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	119 h_u = ops_u.h;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	120 h_v = ops_v.h;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	121
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	122 I_u = speye(m_u);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	123 I_v = speye(m_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	124
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	125 D1_u = ops_u.D1;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	126 H_u = ops_u.H;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	127 Hi_u = ops_u.HI;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	128 e_l_u = ops_u.e_l;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	129 e_r_u = ops_u.e_r;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	130 d1_l_u = ops_u.d1_l;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	131 d1_r_u = ops_u.d1_r;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	132
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	133 D1_v = ops_v.D1;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	134 H_v = ops_v.H;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	135 Hi_v = ops_v.HI;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	136 e_l_v = ops_v.e_l;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	137 e_r_v = ops_v.e_r;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	138 d1_l_v = ops_v.d1_l;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	139 d1_r_v = ops_v.d1_r;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	140
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	141
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	142 % Logical operators
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	143 Du = kr(D1_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	144 Dv = kr(I_u,D1_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	145
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	146 e_w = kr(e_l_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	147 e_e = kr(e_r_u,I_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	148 e_s = kr(I_u,e_l_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	149 e_n = kr(I_u,e_r_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	150
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	151 % Metric coefficients
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	152 coords = g.points();
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	153 x = coords(:,1);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	154 y = coords(:,2);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	155
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	156 x_u = Du*x;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	157 x_v = Dv*x;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	158 y_u = Du*y;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	159 y_v = Dv*y;
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 J = x_u.y_v - x_v.y_u;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	162
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	163 K = cell(dim, dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	164 K{1,1} = y_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	165 K{1,2} = -y_u./J;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	166 K{2,1} = -x_v./J;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	167 K{2,2} = x_u./J;
1306 633757e582e5 Add transformation gradient as property in CurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	168 obj.K = K;
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	169
1211 3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1209 diff changeset	170 % Physical derivatives
3d7faa2ca312 Add physical derivatives Dx and Dy to ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1209 diff changeset	171 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	172 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	173
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	174 % Wrap around Aniosotropic Cartesian
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	175 rho_tilde = J.*rho;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	176
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	177 PHI = cell(dim,dim,dim,dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	178 for i = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	179 for j = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	180 for k = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	181 for l = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	182 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	183 for m = 1:dim
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	184 for n = 1:dim
1209 67eee83fd9c9 Swap indices in anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1205 diff changeset	185 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	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 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	190 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	191 end
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 gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1});
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1295 diff changeset	194 refObj = scheme.Elastic2dVariableAnisotropic(gRef, order, rho_tilde, PHI, opSet, [], hollow);
1205 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 %---- Set object properties ------
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	197 obj.RHO = spdiag(rho);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	198
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	199 % Volume quadrature
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	200 obj.J = spdiag(J);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	201 obj.Ji = spdiag(1./J);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	202 obj.H = obj.J*kr(H_u,H_v);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	203
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	204 % Boundary quadratures
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	205 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	206 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	207 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	208 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	209 obj.s_w = s_w;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	210 obj.s_e = s_e;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	211 obj.s_s = s_s;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	212 obj.s_n = s_n;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	213
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	214 obj.H_w = H_v*spdiag(s_w);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	215 obj.H_e = H_v*spdiag(s_e);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	216 obj.H_s = H_u*spdiag(s_s);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	217 obj.H_n = H_u*spdiag(s_n);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	218
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	219 % Restriction operators
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	220 obj.e_w = refObj.e_w;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	221 obj.e_e = refObj.e_e;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	222 obj.e_s = refObj.e_s;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	223 obj.e_n = refObj.e_n;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	224
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	225 % Adapt things from reference object
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	226 obj.D = refObj.D;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	227 obj.E = refObj.E;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	228
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	229 obj.e1_w = refObj.e1_w;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	230 obj.e1_e = refObj.e1_e;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	231 obj.e1_s = refObj.e1_s;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	232 obj.e1_n = refObj.e1_n;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	233
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	234 obj.e2_w = refObj.e2_w;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	235 obj.e2_e = refObj.e2_e;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	236 obj.e2_s = refObj.e2_s;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	237 obj.e2_n = refObj.e2_n;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	238
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	239 obj.e_scalar_w = refObj.e_scalar_w;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	240 obj.e_scalar_e = refObj.e_scalar_e;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	241 obj.e_scalar_s = refObj.e_scalar_s;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	242 obj.e_scalar_n = refObj.e_scalar_n;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	243
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	244 e1_w = obj.e1_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	245 e1_e = obj.e1_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	246 e1_s = obj.e1_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	247 e1_n = obj.e1_n;
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	248
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	249 e2_w = obj.e2_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	250 e2_e = obj.e2_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	251 e2_s = obj.e2_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	252 e2_n = obj.e2_n;
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	253
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	254 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	255 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	256 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	257 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	258
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	259 obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	260 obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	261 obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	262 obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	263
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	264 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	265 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	266 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	267 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	268
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	269 % Physical normals
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	270 e_w = obj.e_scalar_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	271 e_e = obj.e_scalar_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	272 e_s = obj.e_scalar_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	273 e_n = obj.e_scalar_n;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	274
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	275 e_w_vec = obj.e_w;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	276 e_e_vec = obj.e_e;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	277 e_s_vec = obj.e_s;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	278 e_n_vec = obj.e_n;
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	279
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	280 nu_w = [-1,0];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	281 nu_e = [1,0];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	282 nu_s = [0,-1];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	283 nu_n = [0,1];
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	284
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	285 obj.n_w = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	286 obj.n_e = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	287 obj.n_s = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	288 obj.n_n = cell(2,1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	289
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	290 % 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	291 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	292 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	293 obj.n_w{1} = spdiag(n_w_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	294 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	295 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	296
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	297 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	298 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	299 obj.n_e{1} = spdiag(n_e_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	300 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	301 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	302
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	303 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	304 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	305 obj.n_s{1} = spdiag(n_s_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	306 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	307 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	308
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	309 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	310 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	311 obj.n_n{1} = spdiag(n_n_1);
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	312 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	313 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	314
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	315 % Operators that extract the normal component
43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	316 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	317 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	318 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	319 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	320
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	321 % Operators that extract the tangential component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	322 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	323 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	324 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	325 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	326
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	327 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	328 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	329 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	330 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	331
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	332 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	333 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	334 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	335 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	336
1214 e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	337 % Stress operators
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	338 sigma = cell(dim, dim);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	339 D1 = {obj.Dx, obj.Dy};
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	340 E = obj.E;
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	341 N = length(obj.RHO);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	342 for i = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	343 for j = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	344 sigma{i,j} = sparse(N,2*N);
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	345 for k = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	346 for l = 1:dim
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	347 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	348 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	349 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	350 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	351 end
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	352 obj.sigma = sigma;
e1d4cb8b5309 Add stress operators to scheme anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: 1213 diff changeset	353
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	354 % Misc.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	355 obj.refObj = refObj;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	356 obj.m = refObj.m;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	357 obj.h = refObj.h;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	358 obj.order = order;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	359 obj.grid = g;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	360 obj.dim = dim;
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 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	363
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	364
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	365 % 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	366 % 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	367 % 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	368 % 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	369 % on the first component. Can also be e.g.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	370 % {'normal', 'd'} or {'tangential', 't'} for conditions on
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	371 % tangential/normal component.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	372 % 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	373 % 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	374 % 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	375
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	376 % For displacement bc:
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	377 % bc = {comp, 'd', dComps},
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	378 % where
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	379 % dComps = vector of components with displacement BC. Default: 1:dim.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	380 % 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	381 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	382 default_arg('tuning', 1.0);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	383 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	384
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	385 component = bc{1};
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	386 type = bc{2};
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	387
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	388 switch component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	389
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	390 % If conditions on Cartesian components
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	391 case {1,2}
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	392 [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	393
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	394 % If conditions on normal or tangential components
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	395 case {'n', 't'}
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	396
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	397 switch component
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	398 case 'n'
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	399 en = obj.getBoundaryOperator('en', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	400 case 't'
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	401 en = obj.getBoundaryOperator('et', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	402 end
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	403 e1 = obj.getBoundaryOperator('e1', boundary);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	404
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	405 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	406 [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	407 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	408 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	409
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	410 switch type
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	411 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	412 closure = enen'(c1+c2);
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	413 penalty = ene1'p1;
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	414 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	415 [closure, penalty] = obj.displacementBCNormalTangential(boundary, bc, tuning);
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	416 end
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	417
a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	418 end
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	419
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	420 switch type
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	421 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	422
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	423 s = obj.(['s_' boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	424 s = spdiag(s);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	425 penalty = penalty*s;
1291 a8e730db76e9 Add normal-tangential traction BC ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1272 diff changeset	426
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	427 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	428 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	429
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	430 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	431 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	432 u = obj;
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 component = bc{1};
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	435 type = bc{2};
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	436
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	437 switch component
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	438 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	439 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	440 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	441 N = u.getNormal(boundary);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	442 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	443 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	444 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	445 N = u.getTangent(boundary);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	446 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	447
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	448 % 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	449 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	450 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	451 n = u.getNormal(boundary);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	452
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	453 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	454 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	455 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	456 m_tot = u.grid.N();
1294 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 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	459 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	460
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	461 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	462 dim = u.dim;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	463
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	464 % Preallocate
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	465 [~, 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	466 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	467 penalty = sparse(dim*m_tot, m_int);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	468
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	469 % 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	470 Z = sparse(m_int, m_int);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	471 for i = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	472 for j = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	473 for l = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	474 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	475 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	476 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	477 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	478 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	479 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	480
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	481 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	482 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	483 penalty = penalty - enH_gammaZ;
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	484
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	485 % 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	486 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	487 penalty = penalty - tau_n*H_gamma;
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	488
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	489 % Multiply all terms by inverse of density x quadraure
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	490 closure = RHOiHiclosure;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	491 penalty = RHOiHipenalty;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	492 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	493
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	494 % type Struct that specifies the interface coupling.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	495 % Fields:
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	496 % -- tuning: penalty strength, defaults to 1.0
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	497 % -- interpolation: type of interpolation, default 'none'
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	498 function [closure, penalty, forcingPenalties] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	499
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	500 defaultType.tuning = 1.0;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	501 defaultType.interpolation = 'none';
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	502 defaultType.type = 'standard';
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	503 default_struct('type', defaultType);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	504
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	505 forcingPenalties = [];
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	506
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	507 switch type.type
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	508 case 'standard'
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	509 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	510 case 'normalTangential'
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	511 [closure, penalty, forcingPenalties] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	512 case 'frictionalFault'
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	513 [closure, penalty, forcingPenalties] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	514 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	515
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	516 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	517
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	518 function [closure, penalty, forcingPenalties] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	519 tuning = type.tuning;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	520
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	521 % u denotes the solution in the own domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	522 % v denotes the solution in the neighbour domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	523
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	524 forcingPenalties = cell(1, 1);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	525 u = obj;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	526 v = neighbour_scheme;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	527
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	528 % Operators, u side
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	529 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	530 en_u = u.getBoundaryOperator('en', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	531 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	532 h11_u = u.getBorrowing(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	533 n_u = u.getNormal(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	534
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	535 C_u = u.C;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	536 Ji_u = u.Ji;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	537 s_u = spdiag(u.(['s_' boundary]));
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	538 m_tot_u = u.grid.N();
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	539
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	540 % Operators, v side
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	541 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	542 en_v = v.getBoundaryOperator('en', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	543 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	544 h11_v = v.getBorrowing(neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	545 n_v = v.getNormal(neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	546
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	547 C_v = v.C;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	548 Ji_v = v.Ji;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	549 s_v = spdiag(v.(['s_' neighbour_boundary]));
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	550 m_tot_v = v.grid.N();
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 % Operators that are only required for own domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	553 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	554 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	555
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	556 % Shared operators
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	557 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	558 dim = u.dim;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	559
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	560 % Preallocate
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	561 [~, m_int] = size(H_gamma);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	562 closure = sparse(dimm_tot_u, dimm_tot_u);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	563 penalty = sparse(dimm_tot_u, dimm_tot_v);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	564
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	565 % 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	566 Z_u = sparse(m_int, m_int);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	567 Z_v = sparse(m_int, m_int);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	568 for i = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	569 for j = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	570 for l = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	571 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	572 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	573 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	574 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	575 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	576 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	577 end
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 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	580 closure = closure + en_uH_gammaZ*en_u';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	581 penalty = penalty + en_uH_gammaZ*en_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 % 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	584 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	585 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	586
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	587 % Continuity of normal traction
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	588 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	589 penalty = penalty + 1/2en_uH_gamma*tau_n_v';
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	590 forcing_tau_n = 1/2en_uH_gamma;
1294 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 % 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	593 closure = RHOiHiclosure;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	594 penalty = RHOiHipenalty;
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	595 forcing_tau_n = RHOiHiforcing_tau_n;
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	596
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	597 % ---- Tangential tractions are imposed just like traction BC ------
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	598 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	599
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	600 forcingPenalties{1} = forcing_tau_n;
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	601
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	602 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	603
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	604 % Same interface conditions as in interfaceStandard, but imposed in the normal-tangential
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	605 % coordinate system. For the isotropic case, the components decouple nicely.
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	606 % The resulting scheme is not identical to that of interfaceStandard. This appears to be better.
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	607 function [closure, penalty, forcingPenalties] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type)
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	608 tuning = type.tuning;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	609
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	610 % u denotes the solution in the own domain
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	611 % v denotes the solution in the neighbour domain
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	612
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	613 forcingPenalties = cell(2, 1);
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	614 u = obj;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	615 v = neighbour_scheme;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	616
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	617 % Operators, u side
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	618 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	619 en_u = u.getBoundaryOperator('en', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	620 et_u = u.getBoundaryOperator('et', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	621 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	622 tau_t_u = u.getBoundaryOperator('tau_t', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	623 h11_u = u.getBorrowing(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	624 n_u = u.getNormal(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	625 t_u = u.getTangent(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	626
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	627 C_u = u.C;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	628 Ji_u = u.Ji;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	629 s_u = spdiag(u.(['s_' boundary]));
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	630 m_tot_u = u.grid.N();
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	631
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	632 % Operators, v side
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	633 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	634 en_v = v.getBoundaryOperator('en', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	635 et_v = v.getBoundaryOperator('et', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	636 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	637 tau_t_v = v.getBoundaryOperator('tau_t', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	638 h11_v = v.getBorrowing(neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	639 n_v = v.getNormal(neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	640 t_v = v.getTangent(neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	641
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	642 C_v = v.C;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	643 Ji_v = v.Ji;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	644 s_v = spdiag(v.(['s_' neighbour_boundary]));
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	645 m_tot_v = v.grid.N();
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	646
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	647 % Operators that are only required for own domain
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	648 Hi = u.E{1}inv(u.H)u.E{1}' + u.E{2}inv(u.H)u.E{2}';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	649 RHOi = u.E{1}inv(u.RHO)u.E{1}' + u.E{2}inv(u.RHO)u.E{2}';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	650
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	651 % Shared operators
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	652 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	653 dim = u.dim;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	654
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	655 % Preallocate
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	656 [~, m_int] = size(H_gamma);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	657 closure = sparse(dimm_tot_u, dimm_tot_u);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	658 penalty = sparse(dimm_tot_u, dimm_tot_v);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	659
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	660 % -- Continuity of displacement, term 1: The symmetric term ---
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	661 Zn_u = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	662 Zn_v = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	663 Zt_u = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	664 Zt_v = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	665 for i = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	666 for j = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	667 for l = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	668 for k = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	669 % Penalty strength for normal component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	670 Zn_u = Zn_u + n_u{i}n_u{j}e_u'Ji_uC_u{i,j,k,l}e_un_u{k}*n_u{l};
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	671 Zn_v = Zn_v + n_v{i}n_v{j}e_v'Ji_vC_v{i,j,k,l}e_vn_v{k}*n_v{l};
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	672
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	673 % Penalty strength for tangential component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	674 Zt_u = Zt_u + n_u{i}t_u{j}e_u'Ji_uC_u{i,j,k,l}e_un_u{k}*t_u{l};
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	675 Zt_v = Zt_v + n_v{i}t_v{j}e_v'Ji_vC_v{i,j,k,l}e_vn_v{k}*t_v{l};
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	676 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	677 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	678 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	679 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	680
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	681 Zn = -tuningdim( 1/(4h11_u)s_uZn_u + 1/(4h11_v)s_vZn_v );
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	682 Zt = -tuningdim( 1/(4h11_u)s_uZt_u + 1/(4h11_v)s_vZt_v );
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	683
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	684 % Continuity of normal component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	685 closure = closure + en_uH_gammaZn*en_u';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	686 penalty = penalty + en_uH_gammaZn*en_v';
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	687 forcing_u_n = -en_uH_gammaZn;
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	688
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	689 % Continuity of tangential component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	690 closure = closure + et_uH_gammaZt*et_u';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	691 penalty = penalty + et_uH_gammaZt*et_v';
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	692 forcing_u_t = -et_uH_gammaZt;
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	693 %------------------------------------------------------------------
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	694
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	695 % --- Continuity of displacement, term 2: The symmetrizing term
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	696
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	697 % Continuity of normal displacement
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	698 closure = closure + 1/2tau_n_uH_gamma*en_u';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	699 penalty = penalty + 1/2tau_n_uH_gamma*en_v';
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	700 forcing_u_n = forcing_u_n - 1/2tau_n_uH_gamma;
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	701
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	702 % Continuity of tangential displacement
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	703 closure = closure + 1/2tau_t_uH_gamma*et_u';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	704 penalty = penalty + 1/2tau_t_uH_gamma*et_v';
1319 8b1110385ee2 Bugfix interface forcing in ElasticCurveAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1317 diff changeset	705 forcing_u_t = forcing_u_t - 1/2tau_t_uH_gamma;
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	706 % ------------------------------------------------------------------
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	707
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	708 % --- Continuity of tractions -----------------------------
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	709
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	710 % Continuity of normal traction
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	711 closure = closure - 1/2en_uH_gamma*tau_n_u';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	712 penalty = penalty + 1/2en_uH_gamma*tau_n_v';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	713
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	714 % Continuity of tangential traction
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	715 closure = closure - 1/2et_uH_gamma*tau_t_u';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	716 penalty = penalty + 1/2et_uH_gamma*tau_t_v';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	717 %--------------------------------------------------------------------
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	718
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	719 % Multiply all terms by inverse of density x quadraure
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	720 closure = RHOiHiclosure;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	721 penalty = RHOiHipenalty;
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	722 forcing_u_n = RHOiHiforcing_u_n;
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	723 forcing_u_t = RHOiHiforcing_u_t;
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	724
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	725 forcingPenalties{1} = forcing_u_n;
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	726 forcingPenalties{2} = forcing_u_t;
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	727
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	728 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	729
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	730
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	731 % Returns h11 for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	732 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	733 function h11 = getBorrowing(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	734 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	735
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	736 switch boundary
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	737 case {'w','e'}
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	738 h11 = obj.refObj.h11{1};
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	739 case {'s', 'n'}
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	740 h11 = obj.refObj.h11{2};
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	741 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	742 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	743
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	744 % 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	745 % 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	746 function n = getNormal(obj, boundary)
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	747 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	748
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	749 n = obj.(['n_' boundary]);
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	750 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	751
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	752 % 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	753 % 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	754 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	755 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	756
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	757 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	758 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	759
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	760 % 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	761 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	762 function o = getBoundaryOperator(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	763 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	764 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	765
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	766 o = obj.([op, '_', boundary]);
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	767
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	768 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	769
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	770 % 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	771 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	772 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	773 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	774 assertIsMember(op, {'e'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	775
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	776 switch op
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	777
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	778 case 'e'
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	779 o = obj.(['e_scalar', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	780 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	781
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	782 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	783
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	784 % 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	785 % Formula: tau_i = T_ij u_j
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	786 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	787 function T = getBoundaryTractionOperator(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	788 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	789
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	790 T = obj.(['T', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	791 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	792
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	793 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	794 % corresponding to the number of boundary unknowns
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	795 %
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	796 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	797 function H = getBoundaryQuadrature(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	798 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	799
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	800 H = obj.getBoundaryQuadratureForScalarField(boundary);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	801 I_dim = speye(obj.dim, obj.dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	802 H = kron(H, I_dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	803 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	804
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	805 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	806 % corresponding to the number of boundary grid points
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	807 %
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	808 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	809 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	810 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	811
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	812 H_b = obj.(['H_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	813 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	814
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	815 function N = size(obj)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	816 N = obj.dim*prod(obj.m);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	817 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	818 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	819 end

Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1319:8b1110385ee2 feature/poroelastic