sbplib: +scheme/Elastic2dCurvilinearAnisotropic.m annotate

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1347:ac54767ae1fb feature/poroelastic tip

Add interface, not fully compatible.

author	Martin Almquist <martin.almquist@it.uu.se>
date	Tue, 30 Apr 2024 14:58:35 +0200
parents	412b8ceafbc6
children	df8c71b80c33

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)
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	431 u = obj;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	432
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	433 component = bc{1};
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	434 type = bc{2};
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	435
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	436 switch component
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	437 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	438 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	439 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	440 N = u.getNormal(boundary);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	441 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	442 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	443 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	444 N = u.getTangent(boundary);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	445 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	446
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	447 % 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	448 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	449 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	450 n = u.getNormal(boundary);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	451
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	452 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	453 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	454 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	455 m_tot = u.grid.N();
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	456
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	457 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	458 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	459
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	460 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	461 dim = u.dim;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	462
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	463 % Preallocate
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	464 [~, 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	465 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	466 penalty = sparse(dim*m_tot, m_int);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	467
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	468 % 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	469 Z = sparse(m_int, m_int);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	470 for i = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	471 for j = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	472 for l = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	473 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	474 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	475 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	476 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	477 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
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	480 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	481 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	482 penalty = penalty - enH_gammaZ;
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	483
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	484 % 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	485 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	486 penalty = penalty - tau_n*H_gamma;
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	487
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	488 % 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	489 closure = RHOiHiclosure;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	490 penalty = RHOiHipenalty;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	491 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	492
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	493 % type Struct that specifies the interface coupling.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	494 % Fields:
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	495 % -- tuning: penalty strength, defaults to 1.0
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	496 % -- interpolation: type of interpolation, default 'none'
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	497 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	498
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	499 defaultType.tuning = 1.0;
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	500 defaultType.interpolation = 'none';
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	501 defaultType.type = 'standard';
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	502 default_struct('type', defaultType);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	503
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	504 forcingPenalties = [];
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	505
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	506 switch type.type
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	507 case 'standard'
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	508 [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	509 case 'normalTangential'
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	510 [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	511 case 'frictionalFault'
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	512 [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	513 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	514
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	515 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	516
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	517 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	518 tuning = type.tuning;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	519
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	520 % u denotes the solution in the own domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	521 % v denotes the solution in the neighbour domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	522
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	523 forcingPenalties = cell(1, 1);
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	524 u = obj;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	525 v = neighbour_scheme;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	526
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	527 % Operators, u side
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	528 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	529 en_u = u.getBoundaryOperator('en', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	530 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	531 h11_u = u.getBorrowing(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	532 n_u = u.getNormal(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	533
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	534 C_u = u.C;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	535 Ji_u = u.Ji;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	536 s_u = spdiag(u.(['s_' boundary]));
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	537 m_tot_u = u.grid.N();
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	538
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	539 % Operators, v side
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	540 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	541 en_v = v.getBoundaryOperator('en', neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	542 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	543 h11_v = v.getBorrowing(neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	544 n_v = v.getNormal(neighbour_boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	545
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	546 C_v = v.C;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	547 Ji_v = v.Ji;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	548 s_v = spdiag(v.(['s_' neighbour_boundary]));
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	549 m_tot_v = v.grid.N();
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	550
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	551 % Operators that are only required for own domain
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	552 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	553 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	554
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	555 % Shared operators
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	556 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	557 dim = u.dim;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	558
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	559 % Preallocate
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	560 [~, m_int] = size(H_gamma);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	561 closure = sparse(dimm_tot_u, dimm_tot_u);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	562 penalty = sparse(dimm_tot_u, dimm_tot_v);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	563
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	564 % 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	565 Z_u = sparse(m_int, m_int);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	566 Z_v = sparse(m_int, m_int);
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	567 for i = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	568 for j = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	569 for l = 1:dim
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	570 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	571 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	572 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	573 end
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
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	578 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	579 closure = closure + en_uH_gammaZ*en_u';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	580 penalty = penalty + en_uH_gammaZ*en_v';
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	581
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	582 % 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	583 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	584 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	585
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	586 % Continuity of normal traction
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	587 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	588 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	589 forcing_tau_n = 1/2en_uH_gamma;
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	590
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	591 % 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	592 closure = RHOiHiclosure;
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	593 penalty = RHOiHipenalty;
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	594 forcing_tau_n = RHOiHiforcing_tau_n;
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	595
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	596 % ---- Tangential tractions are imposed just like traction BC ------
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	597 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	598
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	599 forcingPenalties{1} = forcing_tau_n;
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	600
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	601 end
fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	602
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	603 % 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	604 % 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	605 % The resulting scheme is not identical to that of interfaceStandard. This appears to be better.
1322 412b8ceafbc6 Add Zt to output args for Elastic2dCurvilinearAnisotropic.interfaceNormalTangential Martin Almquist <malmquist@stanford.edu> parents: 1321 diff changeset	606 function [closure, penalty, forcingPenalties, Zt] = 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	607 tuning = type.tuning;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	608
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	609 % 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	610 % 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	611
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	612 forcingPenalties = cell(2, 1);
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	613 u = obj;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	614 v = neighbour_scheme;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	615
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	616 % Operators, u side
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	617 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	618 en_u = u.getBoundaryOperator('en', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	619 et_u = u.getBoundaryOperator('et', boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	620 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	621 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	622 h11_u = u.getBorrowing(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	623 n_u = u.getNormal(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	624 t_u = u.getTangent(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	625
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	626 C_u = u.C;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	627 Ji_u = u.Ji;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	628 s_u = spdiag(u.(['s_' boundary]));
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	629 m_tot_u = u.grid.N();
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	630
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	631 % Operators, v side
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	632 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	633 en_v = v.getBoundaryOperator('en', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	634 et_v = v.getBoundaryOperator('et', neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	635 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	636 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	637 h11_v = v.getBorrowing(neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	638 n_v = v.getNormal(neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	639 t_v = v.getTangent(neighbour_boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	640
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	641 C_v = v.C;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	642 Ji_v = v.Ji;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	643 s_v = spdiag(v.(['s_' neighbour_boundary]));
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	644 m_tot_v = v.grid.N();
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	645
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	646 % 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	647 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	648 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	649
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	650 % Shared operators
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	651 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	652 dim = u.dim;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	653
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	654 % Preallocate
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	655 [~, m_int] = size(H_gamma);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	656 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	657 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	658
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	659 % -- 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	660 Zn_u = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	661 Zn_v = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	662 Zt_u = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	663 Zt_v = sparse(m_int, m_int);
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	664 for i = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	665 for j = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	666 for l = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	667 for k = 1:dim
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	668 % Penalty strength for normal component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	669 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	670 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	671
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	672 % Penalty strength for tangential component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	673 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	674 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	675 end
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
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	680 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	681 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	682
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	683 % Continuity of normal component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	684 closure = closure + en_uH_gammaZn*en_u';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	685 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	686 forcing_u_n = -en_uH_gammaZn;
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	687
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	688 % Continuity of tangential component
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	689 closure = closure + et_uH_gammaZt*et_u';
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	690 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	691 forcing_u_t = -et_uH_gammaZt;
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	692 %------------------------------------------------------------------
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 % --- 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	695
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	696 % Continuity of normal displacement
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	697 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	698 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	699 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	700
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	701 % Continuity of tangential displacement
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	702 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	703 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	704 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	705 % ------------------------------------------------------------------
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 % --- Continuity of tractions -----------------------------
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	708
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	709 % Continuity of normal traction
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	710 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	711 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	712
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	713 % Continuity of tangential traction
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	714 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	715 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	716 %--------------------------------------------------------------------
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 % 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	719 closure = RHOiHiclosure;
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	720 penalty = RHOiHipenalty;
1317 34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	721 forcing_u_n = RHOiHiforcing_u_n;
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	722 forcing_u_t = RHOiHiforcing_u_t;
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	723
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	724 forcingPenalties{1} = forcing_u_n;
34997aced843 Add some interface forcing penalties in ElasticCurvilinearAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1315 diff changeset	725 forcingPenalties{2} = forcing_u_t;
1315 6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	726
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	727 end
6c308da9dcbc Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic. Martin Almquist <malmquist@stanford.edu> parents: 1306 diff changeset	728
1294 fe712a1fca3f Add frictional fault interface in Elastic2dAnisotropicCurve Martin Almquist <malmquist@stanford.edu> parents: 1291 diff changeset	729
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	730 % Returns h11 for the boundary specified by the string boundary.
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	731 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	732 function h11 = getBorrowing(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	733 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	734
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	735 switch boundary
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	736 case {'w','e'}
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	737 h11 = obj.refObj.h11{1};
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	738 case {'s', 'n'}
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	739 h11 = obj.refObj.h11{2};
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	740 end
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
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	743 % 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	744 % 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	745 function n = getNormal(obj, boundary)
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	746 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	747
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	748 n = obj.(['n_' boundary]);
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	749 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	750
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	751 % 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	752 % 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	753 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	754 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	755
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 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	757 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	758
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	759 % 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	760 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	761 function o = getBoundaryOperator(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	762 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	763 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	764
1213 43f1cd11e8e8 Add physical normals to AnisotropicCurvilinear Martin Almquist <malmquist@stanford.edu> parents: 1211 diff changeset	765 o = obj.([op, '_', boundary]);
1205 3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	766
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	767 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	768
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	769 % 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	770 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	771 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	772 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	773 assertIsMember(op, {'e'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	774
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	775 switch op
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	776
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	777 case 'e'
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	778 o = obj.(['e_scalar', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	779 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	780
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	781 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	782
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	783 % 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	784 % Formula: tau_i = T_ij u_j
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	785 % op -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	786 function T = getBoundaryTractionOperator(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	787 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	788
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	789 T = obj.(['T', '_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	790 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	791
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	792 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	793 % corresponding to the number of boundary unknowns
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	794 %
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	795 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	796 function H = getBoundaryQuadrature(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	797 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	798
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	799 H = obj.getBoundaryQuadratureForScalarField(boundary);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	800 I_dim = speye(obj.dim, obj.dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	801 H = kron(H, I_dim);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	802 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	803
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	804 % Returns square boundary quadrature matrix, of dimension
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	805 % corresponding to the number of boundary grid points
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	806 %
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	807 % boundary -- string
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	808 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	809 assertIsMember(boundary, {'w', 'e', 's', 'n'})
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	810
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	811 H_b = obj.(['H_', boundary]);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	812 end
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	813
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	814 function N = size(obj)
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	815 N = obj.dim*prod(obj.m);
3258dca12af8 Add scheme for anisotropic curvilinear Martin Almquist <malmquist@stanford.edu> parents: diff changeset	816 end
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

Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1347:ac54767ae1fb feature/poroelastic tip