sbplib: +scheme/Elastic2dCurvilinearAnisotropic.m annotate

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1294:fe712a1fca3f feature/poroelastic

Add frictional fault interface in Elastic2dAnisotropicCurve

author	Martin Almquist <malmquist@stanford.edu>
date	Thu, 02 Jul 2020 15:00:21 -0700
parents	a8e730db76e9
children	cb053fabbedc

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

Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinearAnisotropic.m @ 1294:fe712a1fca3f feature/poroelastic