sbplib: +scheme/Elastic2dCurvilinear.m annotate

author	Martin Almquist <malmquist@stanford.edu>
date	Mon, 07 May 2018 14:35:54 -0700
parents
children	f4e2a6a2df08

rev	line source
739 8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	1 classdef Elastic2dCurvilinear < scheme.Scheme
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	2
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	3 % Discretizes the elastic wave equation in curvilinear coordinates.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	4 %
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	5 % Untransformed equation:
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	6 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	7 %
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	8 % Transformed equation:
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	9 % J*rho u_{i,tt} = dk J b_ik lambda b_jl dl u_j
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	10 % + dk J b_jk mu b_il dl u_j
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	11 % + dk J b_jk mu b_jl dl u_i
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	12 % opSet should be cell array of opSets, one per dimension. This
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	13 % is useful if we have periodic BC in one direction.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	14
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	15 properties
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	16 m % Number of points in each direction, possibly a vector
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	17 h % Grid spacing
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	18
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	19 grid
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	20 dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	21
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	22 order % Order of accuracy for the approximation
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	23
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	24 % Diagonal matrices for varible coefficients
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	25 LAMBDA % Variable coefficient, related to dilation
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	26 MU % Shear modulus, variable coefficient
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	27 RHO, RHOi % Density, variable
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	28
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	29 % Metric coefficients
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	30 b % Cell matrix of size dim x dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	31 J, Ji
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	32
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	33 D % Total operator
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	34 D1 % First derivatives
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	35
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	36 % Second derivatives
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	37 D2_lambda
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	38 D2_mu
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	39
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	40 % Traction operators used for BC
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	41 T_l, T_r
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	42 tau_l, tau_r
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	43
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	44 H, Hi % Inner products
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	45 phi % Borrowing constant for (d1 - e^T*D1) from R
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	46 gamma % Borrowing constant for d1 from M
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	47 H11 % First element of H
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	48 e_l, e_r
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	49 d1_l, d1_r % Normal derivatives at the boundary
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	50 E % E{i}^T picks out component i
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	51
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	52 H_boundary_l, H_boundary_r % Boundary inner products
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	53
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	54 % Kroneckered norms and coefficients
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	55 RHOi_kron
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	56 Ji_kron, J_kron
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	57 Hi_kron, H_kron
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	58 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	59
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	60 methods
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	61
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	62 function obj = Elastic2dCurvilinear(g ,order, lambda_fun, mu_fun, rho_fun, opSet)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	63 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	64 default_arg('lambda_fun', @(x,y) 0*x+1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	65 default_arg('mu_fun', @(x,y) 0*x+1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	66 default_arg('rho_fun', @(x,y) 0*x+1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	67 dim = 2;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	68
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	69 lambda = grid.evalOn(g, lambda_fun);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	70 mu = grid.evalOn(g, mu_fun);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	71 rho = grid.evalOn(g, rho_fun);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	72 m = g.size();
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	73 obj.m = m;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	74 m_tot = g.N();
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	75
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	76 % 1D operators
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	77 ops = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	78 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	79 ops{i} = opSet{i}(m(i), {0, 1}, order);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	80 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	81
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	82 % Borrowing constants
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	83 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	84 beta = ops{i}.borrowing.R.delta_D;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	85 obj.H11{i} = ops{i}.borrowing.H11;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	86 obj.phi{i} = beta/obj.H11{i};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	87 obj.gamma{i} = ops{i}.borrowing.M.d1;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	88 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	89
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	90 I = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	91 D1 = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	92 D2 = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	93 H = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	94 Hi = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	95 e_l = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	96 e_r = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	97 d1_l = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	98 d1_r = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	99
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	100 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	101 I{i} = speye(m(i));
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	102 D1{i} = ops{i}.D1;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	103 D2{i} = ops{i}.D2;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	104 H{i} = ops{i}.H;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	105 Hi{i} = ops{i}.HI;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	106 e_l{i} = ops{i}.e_l;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	107 e_r{i} = ops{i}.e_r;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	108 d1_l{i} = ops{i}.d1_l;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	109 d1_r{i} = ops{i}.d1_r;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	110 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	111
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	112 %====== Assemble full operators ========
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	113
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	114 % Variable coefficients
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	115 LAMBDA = spdiag(lambda);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	116 obj.LAMBDA = LAMBDA;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	117 MU = spdiag(mu);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	118 obj.MU = MU;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	119 RHO = spdiag(rho);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	120 obj.RHO = RHO;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	121 obj.RHOi = inv(RHO);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	122
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	123 % Allocate
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	124 obj.D1 = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	125 obj.D2_lambda = cell(dim,dim,dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	126 obj.D2_mu = cell(dim,dim,dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	127 obj.e_l = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	128 obj.e_r = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	129 obj.d1_l = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	130 obj.d1_r = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	131
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	132 % D1
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	133 obj.D1{1} = kron(D1{1},I{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	134 obj.D1{2} = kron(I{1},D1{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	135
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	136 % -- Metric coefficients ----
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	137 coords = g.points();
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	138 x = coords(:,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	139 y = coords(:,2);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	140
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	141 x_xi = obj.D1{1}*x;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	142 x_eta = obj.D1{2}*x;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	143 y_xi = obj.D1{1}*y;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	144 y_eta = obj.D1{2}*y;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	145
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	146 J = x_xi.y_eta - x_eta.y_xi;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	147
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	148 b = cell(dim,dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	149 b{1,1} = y_eta./J;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	150 b{1,2} = -x_eta./J;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	151 b{2,1} = -y_xi./J;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	152 b{2,2} = x_xi./J;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	153
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	154 % Scale factors for boundary integrals
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	155 beta = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	156 beta{1} = sqrt(x_eta.^2 + y_eta.^2);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	157 beta{2} = sqrt(x_xi.^2 + y_xi.^2);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	158
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	159
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	160 J = spdiag(J);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	161 Ji = inv(J);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	162 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	163 beta{i} = spdiag(beta{i});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	164 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	165 b{i,j} = spdiag(b{i,j});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	166 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	167 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	168 obj.J = J;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	169 obj.Ji = Ji;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	170 obj.b = b;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	171 %----------------------------
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	172
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	173 % Boundary operators
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	174 obj.e_l{1} = kron(e_l{1},I{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	175 obj.e_l{2} = kron(I{1},e_l{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	176 obj.e_r{1} = kron(e_r{1},I{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	177 obj.e_r{2} = kron(I{1},e_r{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	178
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	179 obj.d1_l{1} = kron(d1_l{1},I{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	180 obj.d1_l{2} = kron(I{1},d1_l{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	181 obj.d1_r{1} = kron(d1_r{1},I{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	182 obj.d1_r{2} = kron(I{1},d1_r{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	183
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	184 % D2
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	185 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	186 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	187 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	188 obj.D2_lambda{i,j,k} = sparse(m_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	189 obj.D2_mu{i,j,k} = sparse(m_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	190 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	191 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	192 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	193 ind = grid.funcToMatrix(g, 1:m_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	194
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	195 % x-dir
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	196 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	197 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	198 for k = 1
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	199
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	200 coeff_lambda = Jb{i,k}b{j,k}*lambda;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	201 coeff_mu = Jb{j,k}b{i,k}*mu;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	202
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	203 for col = 1:m(2)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	204 D_lambda = D2{1}(coeff_lambda(ind(:,col)));
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	205 D_mu = D2{1}(coeff_mu(ind(:,col)));
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	206
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	207 p = ind(:,col);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	208 obj.D2_lambda{i,j,k}(p,p) = D_lambda;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	209 obj.D2_mu{i,j,k}(p,p) = D_mu;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	210 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	211
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	212 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	213 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	214 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	215
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	216 % y-dir
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	217 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	218 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	219 for k = 2
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	220
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	221 coeff_lambda = Jb{i,k}b{j,k}*lambda;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	222 coeff_mu = Jb{j,k}b{i,k}*mu;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	223
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	224 for row = 1:m(1)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	225 D_lambda = D2{2}(coeff_lambda(ind(row,:)));
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	226 D_mu = D2{2}(coeff_mu(ind(row,:)));
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	227
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	228 p = ind(row,:);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	229 obj.D2_lambda{i,j,k}(p,p) = D_lambda;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	230 obj.D2_mu{i,j,k}(p,p) = D_mu;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	231 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	232
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	233 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	234 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	235 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	236
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	237 % Quadratures
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	238 obj.H = kron(H{1},H{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	239 obj.Hi = inv(obj.H);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	240 obj.H_boundary_l = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	241 obj.H_boundary_l{1} = obj.e_l{1}'beta{1}obj.e_l{1}*H{2};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	242 obj.H_boundary_l{2} = obj.e_l{2}'beta{2}obj.e_l{2}*H{1};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	243 obj.H_boundary_r = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	244 obj.H_boundary_r{1} = obj.e_r{1}'beta{1}obj.e_r{1}*H{2};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	245 obj.H_boundary_r{2} = obj.e_r{2}'beta{2}obj.e_r{2}*H{1};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	246
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	247 % E{i}^T picks out component i.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	248 E = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	249 I = speye(m_tot,m_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	250 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	251 e = sparse(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	252 e(i) = 1;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	253 E{i} = kron(I,e);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	254 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	255 obj.E = E;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	256
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	257 % Differentiation matrix D (without SAT)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	258 D2_lambda = obj.D2_lambda;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	259 D2_mu = obj.D2_mu;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	260 D1 = obj.D1;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	261 D = sparse(dimm_tot,dimm_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	262 d = @kroneckerDelta; % Kronecker delta
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	263 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	264 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	265 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	266 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	267 for l = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	268 D = D + E{i}Jiinv(RHO)( d(k,l)D2_lambda{i,j,k}*E{j}' + ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	269 db(k,l)D1{k}Jb{i,k}b{j,l}LAMBDAD1{l}*E{j}' ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	270 );
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	271
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	272 D = D + E{i}Jiinv(RHO)( d(k,l)D2_mu{i,j,k}*E{j}' + ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	273 db(k,l)D1{k}Jb{j,k}b{i,l}MUD1{l}*E{j}' ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	274 );
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	275
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	276 D = D + E{i}Jiinv(RHO)( d(k,l)D2_mu{j,j,k}*E{i}' + ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	277 db(k,l)D1{k}Jb{j,k}b{j,l}MUD1{l}*E{i}' ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	278 );
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	279
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	280 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	281 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	282 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	283 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	284 obj.D = D;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	285 %=========================================%
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	286
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	287 % Numerical traction operators for BC.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	288 % Because d1 =/= e0^T*D1, the numerical tractions are different
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	289 % at every boundary.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	290 T_l = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	291 T_r = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	292 tau_l = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	293 tau_r = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	294 % tau^{j}_i = sum_k T^{j}_{ik} u_k
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	295
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	296 d1_l = obj.d1_l;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	297 d1_r = obj.d1_r;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	298 e_l = obj.e_l;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	299 e_r = obj.e_r;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	300
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	301 % Loop over boundaries
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	302 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	303 T_l{j} = cell(dim,dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	304 T_r{j} = cell(dim,dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	305 tau_l{j} = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	306 tau_r{j} = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	307
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	308 % Loop over components
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	309 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	310 tau_l{j}{i} = sparse(m_tot,dim*m_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	311 tau_r{j}{i} = sparse(m_tot,dim*m_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	312
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	313 % Loop over components that T_{ik}^{(j)} acts on
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	314 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	315
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	316 T_l{j}{i,k} = sparse(m_tot,m_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	317 T_r{j}{i,k} = sparse(m_tot,m_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	318
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	319 for m = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	320 for l = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	321 T_l{j}{i,k} = T_l{j}{i,k} + ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	322 -d(k,l)* Jb{i,j}b{k,m}LAMBDA(d(m,j)e_l{m}d1_l{m}' + db(m,j)*D1{m}) ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	323 -d(k,l)* Jb{k,j}b{i,m}MU(d(m,j)e_l{m}d1_l{m}' + db(m,j)*D1{m}) ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	324 -d(i,k)* Jb{l,j}b{l,m}MU(d(m,j)e_l{m}d1_l{m}' + db(m,j)*D1{m});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	325
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	326 T_r{j}{i,k} = T_r{j}{i,k} + ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	327 d(k,l)* Jb{i,j}b{k,m}LAMBDA(d(m,j)e_r{m}d1_r{m}' + db(m,j)*D1{m}) + ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	328 d(k,l)* Jb{k,j}b{i,m}MU(d(m,j)e_r{m}d1_r{m}' + db(m,j)*D1{m}) + ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	329 d(i,k)* Jb{l,j}b{l,m}MU(d(m,j)e_r{m}d1_r{m}' + db(m,j)*D1{m});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	330 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	331 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	332
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	333 T_l{j}{i,k} = inv(beta{j})*T_l{j}{i,k};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	334 T_r{j}{i,k} = inv(beta{j})*T_r{j}{i,k};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	335
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	336 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	337 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	338 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	339
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	340 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	341 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	342 obj.T_l = T_l;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	343 obj.T_r = T_r;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	344 obj.tau_l = tau_l;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	345 obj.tau_r = tau_r;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	346
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	347 % Kroneckered norms and coefficients
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	348 I_dim = speye(dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	349 obj.RHOi_kron = kron(obj.RHOi, I_dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	350 obj.Ji_kron = kron(obj.Ji, I_dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	351 obj.Hi_kron = kron(obj.Hi, I_dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	352 obj.H_kron = kron(obj.H, I_dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	353 obj.J_kron = kron(obj.J, I_dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	354
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	355 % Misc.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	356 obj.h = g.scaling();
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	357 obj.order = order;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	358 obj.grid = g;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	359 obj.dim = dim;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	360
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	361 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	362
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	363
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	364 % Closure functions return the operators applied to the own domain to close the boundary
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	365 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	366 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	367 % type is a cell array of strings specifying the type of boundary condition for each component.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	368 % data is a function returning the data that should be applied at the boundary.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	369 % neighbour_scheme is an instance of Scheme that should be interfaced to.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	370 % neighbour_boundary is a string specifying which boundary to interface to.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	371 function [closure, penalty] = boundary_condition(obj, boundary, type, tuning)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	372 default_arg('type',{'free','free'});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	373 default_arg('tuning', 1.2);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	374
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	375 if ~iscell(type)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	376 type = {type, type};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	377 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	378
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	379 % j is the coordinate direction of the boundary
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	380 j = obj.get_boundary_number(boundary);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	381 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	382
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	383 E = obj.E;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	384 Hi = obj.Hi;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	385 LAMBDA = obj.LAMBDA;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	386 MU = obj.MU;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	387 RHOi = obj.RHOi;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	388 Ji = obj.Ji;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	389
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	390 dim = obj.dim;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	391 m_tot = obj.grid.N();
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	392
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	393 % Preallocate
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	394 closure = sparse(dimm_tot, dimm_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	395 penalty = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	396 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	397 penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j));
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	398 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	399
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	400 % Loop over components that we (potentially) have different BC on
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	401 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	402 switch type{k}
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	403
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	404 % Dirichlet boundary condition
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	405 case {'D','d','dirichlet','Dirichlet'}
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	406
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	407 phi = obj.phi{j};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	408 h = obj.h(j);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	409 h11 = obj.H11{j}*h;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	410 gamma = obj.gamma{j};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	411
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	412 a_lambda = dim/h11 + 1/(h11*phi);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	413 a_mu_i = 2/(gamma*h);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	414 a_mu_ij = 2/h11 + 1/(h11*phi);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	415
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	416 d = @kroneckerDelta; % Kronecker delta
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	417 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	418 alpha = @(i,j) tuning( d(i,j) a_lambda*LAMBDA ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	419 + d(i,j)* a_mu_i*MU ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	420 + db(i,j)a_mu_ijMU );
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	421
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	422 % Loop over components that Dirichlet penalties end up on
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	423 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	424 C = T{k,i};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	425 A = -d(i,k)*alpha(i,j);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	426 B = A + C;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	427 closure = closure + E{i}RHOiHiJiB'eH_gamma(e'E{k}' );
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	428 penalty{k} = penalty{k} - E{i}RHOiHiJiB'eH_gamma;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	429 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	430
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	431 % Free boundary condition
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	432 case {'F','f','Free','free','traction','Traction','t','T'}
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	433 closure = closure - E{k}RHOiJiHieH_gamma (e'*tau{k} );
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	434 penalty{k} = penalty{k} + E{k}RHOiJiHie*H_gamma;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	435
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	436 % Unknown boundary condition
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	437 otherwise
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	438 error('No such boundary condition: type = %s',type);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	439 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	440 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	441 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	442
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	443 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	444 % u denotes the solution in the own domain
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	445 % v denotes the solution in the neighbour domain
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	446 % Operators without subscripts are from the own domain.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	447 error('Not implemented');
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	448 tuning = 1.2;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	449
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	450 % j is the coordinate direction of the boundary
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	451 j = obj.get_boundary_number(boundary);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	452 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	453
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	454 % Get boundary operators
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	455 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	456 [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	457
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	458 % Operators and quantities that correspond to the own domain only
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	459 Hi = obj.Hi;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	460 RHOi = obj.RHOi;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	461 dim = obj.dim;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	462
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	463 %--- Other operators ----
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	464 m_tot_u = obj.grid.N();
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	465 E = obj.E;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	466 LAMBDA_u = obj.LAMBDA;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	467 MU_u = obj.MU;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	468 lambda_u = e'LAMBDA_ue;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	469 mu_u = e'MU_ue;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	470
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	471 m_tot_v = neighbour_scheme.grid.N();
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	472 E_v = neighbour_scheme.E;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	473 LAMBDA_v = neighbour_scheme.LAMBDA;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	474 MU_v = neighbour_scheme.MU;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	475 lambda_v = e_v'LAMBDA_ve_v;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	476 mu_v = e_v'MU_ve_v;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	477 %-------------------------
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	478
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	479 % Borrowing constants
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	480 phi_u = obj.phi{j};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	481 h_u = obj.h(j);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	482 h11_u = obj.H11{j}*h_u;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	483 gamma_u = obj.gamma{j};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	484
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	485 phi_v = neighbour_scheme.phi{j_v};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	486 h_v = neighbour_scheme.h(j_v);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	487 h11_v = neighbour_scheme.H11{j_v}*h_v;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	488 gamma_v = neighbour_scheme.gamma{j_v};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	489
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	490 % E > sum_i 1/(2alpha_ij)(tau_i)^2
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	491 function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	492 th1 = h11/(2*dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	493 th2 = h11*phi/2;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	494 th3 = h*gamma;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	495 a1 = ( (th1 + th2)th3lambda + 4th1th2mu ) / (2th1th2th3);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	496 a2 = ( 16(th1 + th2)lambdamu ) / (th1th2*th3);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	497 alpha_ii = a1 + sqrt(a2 + a1^2);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	498
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	499 alpha_ij = mu(2/h11 + 1/(phih11));
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	500 end

739

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

1 classdef Elastic2dCurvilinear < scheme.Scheme

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

2

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

3 % Discretizes the elastic wave equation in curvilinear coordinates.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

4 %

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

5 % Untransformed equation:

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

6 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

7 %

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

8 % Transformed equation:

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

9 % J*rho u_{i,tt} = dk J b_ik lambda b_jl dl u_j

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

10 % + dk J b_jk mu b_il dl u_j

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

11 % + dk J b_jk mu b_jl dl u_i

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

12 % opSet should be cell array of opSets, one per dimension. This

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

13 % is useful if we have periodic BC in one direction.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

14

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

15 properties

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

16 m % Number of points in each direction, possibly a vector

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

17 h % Grid spacing

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

18

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

19 grid

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

20 dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

21

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

22 order % Order of accuracy for the approximation

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

23

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

24 % Diagonal matrices for varible coefficients

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

25 LAMBDA % Variable coefficient, related to dilation

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

26 MU % Shear modulus, variable coefficient

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

27 RHO, RHOi % Density, variable

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

28

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

29 % Metric coefficients

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

30 b % Cell matrix of size dim x dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

31 J, Ji

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

32

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

33 D % Total operator

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

34 D1 % First derivatives

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

35

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

36 % Second derivatives

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

37 D2_lambda

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

38 D2_mu

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

39

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

40 % Traction operators used for BC

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

41 T_l, T_r

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

42 tau_l, tau_r

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

43

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

44 H, Hi % Inner products

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

45 phi % Borrowing constant for (d1 - e^T*D1) from R

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

46 gamma % Borrowing constant for d1 from M

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

47 H11 % First element of H

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

48 e_l, e_r

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

49 d1_l, d1_r % Normal derivatives at the boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

50 E % E{i}^T picks out component i

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

51

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

52 H_boundary_l, H_boundary_r % Boundary inner products

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

53

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

54 % Kroneckered norms and coefficients

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

55 RHOi_kron

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

56 Ji_kron, J_kron

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

57 Hi_kron, H_kron

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

58 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

59

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

60 methods

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

61

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

62 function obj = Elastic2dCurvilinear(g ,order, lambda_fun, mu_fun, rho_fun, opSet)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

63 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

64 default_arg('lambda_fun', @(x,y) 0*x+1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

65 default_arg('mu_fun', @(x,y) 0*x+1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

66 default_arg('rho_fun', @(x,y) 0*x+1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

67 dim = 2;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

68

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

69 lambda = grid.evalOn(g, lambda_fun);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

70 mu = grid.evalOn(g, mu_fun);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

71 rho = grid.evalOn(g, rho_fun);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

72 m = g.size();

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

73 obj.m = m;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

74 m_tot = g.N();

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

75

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

76 % 1D operators

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

77 ops = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

78 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

79 ops{i} = opSet{i}(m(i), {0, 1}, order);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

80 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

81

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

82 % Borrowing constants

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

83 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

84 beta = ops{i}.borrowing.R.delta_D;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

85 obj.H11{i} = ops{i}.borrowing.H11;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

86 obj.phi{i} = beta/obj.H11{i};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

87 obj.gamma{i} = ops{i}.borrowing.M.d1;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

88 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

89

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

90 I = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

91 D1 = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

92 D2 = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

93 H = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

94 Hi = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

95 e_l = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

96 e_r = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

97 d1_l = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

98 d1_r = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

99

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

100 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

101 I{i} = speye(m(i));

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

102 D1{i} = ops{i}.D1;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

103 D2{i} = ops{i}.D2;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

104 H{i} = ops{i}.H;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

105 Hi{i} = ops{i}.HI;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

106 e_l{i} = ops{i}.e_l;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

107 e_r{i} = ops{i}.e_r;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

108 d1_l{i} = ops{i}.d1_l;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

109 d1_r{i} = ops{i}.d1_r;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

110 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

111

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

112 %====== Assemble full operators ========

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

113

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

114 % Variable coefficients

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

115 LAMBDA = spdiag(lambda);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

116 obj.LAMBDA = LAMBDA;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

117 MU = spdiag(mu);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

118 obj.MU = MU;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

119 RHO = spdiag(rho);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

120 obj.RHO = RHO;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

121 obj.RHOi = inv(RHO);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

122

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

123 % Allocate

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

124 obj.D1 = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

125 obj.D2_lambda = cell(dim,dim,dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

126 obj.D2_mu = cell(dim,dim,dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

127 obj.e_l = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

128 obj.e_r = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

129 obj.d1_l = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

130 obj.d1_r = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

131

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

132 % D1

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

133 obj.D1{1} = kron(D1{1},I{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

134 obj.D1{2} = kron(I{1},D1{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

135

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

136 % -- Metric coefficients ----

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

137 coords = g.points();

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

138 x = coords(:,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

139 y = coords(:,2);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

140

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

141 x_xi = obj.D1{1}*x;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

142 x_eta = obj.D1{2}*x;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

143 y_xi = obj.D1{1}*y;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

144 y_eta = obj.D1{2}*y;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

145

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

146 J = x_xi.*y_eta - x_eta.*y_xi;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

147

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

148 b = cell(dim,dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

149 b{1,1} = y_eta./J;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

150 b{1,2} = -x_eta./J;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

151 b{2,1} = -y_xi./J;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

152 b{2,2} = x_xi./J;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

153

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

154 % Scale factors for boundary integrals

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

155 beta = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

156 beta{1} = sqrt(x_eta.^2 + y_eta.^2);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

157 beta{2} = sqrt(x_xi.^2 + y_xi.^2);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

158

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

159

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

160 J = spdiag(J);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

161 Ji = inv(J);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

162 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

163 beta{i} = spdiag(beta{i});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

164 for j = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

165 b{i,j} = spdiag(b{i,j});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

166 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

167 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

168 obj.J = J;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

169 obj.Ji = Ji;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

170 obj.b = b;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

171 %----------------------------

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

172

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

173 % Boundary operators

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

174 obj.e_l{1} = kron(e_l{1},I{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

175 obj.e_l{2} = kron(I{1},e_l{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

176 obj.e_r{1} = kron(e_r{1},I{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

177 obj.e_r{2} = kron(I{1},e_r{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

178

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

179 obj.d1_l{1} = kron(d1_l{1},I{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

180 obj.d1_l{2} = kron(I{1},d1_l{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

181 obj.d1_r{1} = kron(d1_r{1},I{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

182 obj.d1_r{2} = kron(I{1},d1_r{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

183

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

184 % D2

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

185 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

186 for j = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

187 for k = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

188 obj.D2_lambda{i,j,k} = sparse(m_tot);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

189 obj.D2_mu{i,j,k} = sparse(m_tot);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

190 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

191 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

192 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

193 ind = grid.funcToMatrix(g, 1:m_tot);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

194

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

195 % x-dir

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

196 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

197 for j = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

198 for k = 1

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

199

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

200 coeff_lambda = J*b{i,k}*b{j,k}*lambda;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

201 coeff_mu = J*b{j,k}*b{i,k}*mu;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

202

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

203 for col = 1:m(2)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

204 D_lambda = D2{1}(coeff_lambda(ind(:,col)));

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

205 D_mu = D2{1}(coeff_mu(ind(:,col)));

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

206

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

207 p = ind(:,col);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

208 obj.D2_lambda{i,j,k}(p,p) = D_lambda;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

209 obj.D2_mu{i,j,k}(p,p) = D_mu;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

210 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

211

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

212 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

213 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

214 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

215

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

216 % y-dir

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

217 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

218 for j = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

219 for k = 2

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

220

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

221 coeff_lambda = J*b{i,k}*b{j,k}*lambda;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

222 coeff_mu = J*b{j,k}*b{i,k}*mu;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

223

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

224 for row = 1:m(1)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

225 D_lambda = D2{2}(coeff_lambda(ind(row,:)));

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

226 D_mu = D2{2}(coeff_mu(ind(row,:)));

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

227

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

228 p = ind(row,:);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

229 obj.D2_lambda{i,j,k}(p,p) = D_lambda;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

230 obj.D2_mu{i,j,k}(p,p) = D_mu;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

231 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

232

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

233 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

234 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

235 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

236

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

237 % Quadratures

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

238 obj.H = kron(H{1},H{2});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

239 obj.Hi = inv(obj.H);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

240 obj.H_boundary_l = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

241 obj.H_boundary_l{1} = obj.e_l{1}'*beta{1}*obj.e_l{1}*H{2};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

242 obj.H_boundary_l{2} = obj.e_l{2}'*beta{2}*obj.e_l{2}*H{1};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

243 obj.H_boundary_r = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

244 obj.H_boundary_r{1} = obj.e_r{1}'*beta{1}*obj.e_r{1}*H{2};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

245 obj.H_boundary_r{2} = obj.e_r{2}'*beta{2}*obj.e_r{2}*H{1};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

246

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

247 % E{i}^T picks out component i.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

248 E = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

249 I = speye(m_tot,m_tot);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

250 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

251 e = sparse(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

252 e(i) = 1;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

253 E{i} = kron(I,e);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

254 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

255 obj.E = E;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

256

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

257 % Differentiation matrix D (without SAT)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

258 D2_lambda = obj.D2_lambda;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

259 D2_mu = obj.D2_mu;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

260 D1 = obj.D1;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

261 D = sparse(dim*m_tot,dim*m_tot);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

262 d = @kroneckerDelta; % Kronecker delta

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

263 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

264 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

265 for j = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

266 for k = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

267 for l = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

268 D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_lambda{i,j,k}*E{j}' + ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

269 db(k,l)*D1{k}*J*b{i,k}*b{j,l}*LAMBDA*D1{l}*E{j}' ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

270 );

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

271

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

272 D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_mu{i,j,k}*E{j}' + ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

273 db(k,l)*D1{k}*J*b{j,k}*b{i,l}*MU*D1{l}*E{j}' ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

274 );

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

275

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

276 D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_mu{j,j,k}*E{i}' + ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

277 db(k,l)*D1{k}*J*b{j,k}*b{j,l}*MU*D1{l}*E{i}' ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

278 );

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

279

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

280 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

281 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

282 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

283 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

284 obj.D = D;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

285 %=========================================%

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

286

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

287 % Numerical traction operators for BC.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

288 % Because d1 =/= e0^T*D1, the numerical tractions are different

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

289 % at every boundary.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

290 T_l = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

291 T_r = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

292 tau_l = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

293 tau_r = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

294 % tau^{j}_i = sum_k T^{j}_{ik} u_k

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

295

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

296 d1_l = obj.d1_l;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

297 d1_r = obj.d1_r;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

298 e_l = obj.e_l;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

299 e_r = obj.e_r;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

300

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

301 % Loop over boundaries

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

302 for j = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

303 T_l{j} = cell(dim,dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

304 T_r{j} = cell(dim,dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

305 tau_l{j} = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

306 tau_r{j} = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

307

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

308 % Loop over components

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

309 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

310 tau_l{j}{i} = sparse(m_tot,dim*m_tot);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

311 tau_r{j}{i} = sparse(m_tot,dim*m_tot);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

312

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

313 % Loop over components that T_{ik}^{(j)} acts on

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

314 for k = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

315

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

316 T_l{j}{i,k} = sparse(m_tot,m_tot);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

317 T_r{j}{i,k} = sparse(m_tot,m_tot);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

318

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

319 for m = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

320 for l = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

321 T_l{j}{i,k} = T_l{j}{i,k} + ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

322 -d(k,l)* J*b{i,j}*b{k,m}*LAMBDA*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}) ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

323 -d(k,l)* J*b{k,j}*b{i,m}*MU*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}) ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

324 -d(i,k)* J*b{l,j}*b{l,m}*MU*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

325

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

326 T_r{j}{i,k} = T_r{j}{i,k} + ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

327 d(k,l)* J*b{i,j}*b{k,m}*LAMBDA*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}) + ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

328 d(k,l)* J*b{k,j}*b{i,m}*MU*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}) + ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

329 d(i,k)* J*b{l,j}*b{l,m}*MU*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

330 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

331 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

332

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

333 T_l{j}{i,k} = inv(beta{j})*T_l{j}{i,k};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

334 T_r{j}{i,k} = inv(beta{j})*T_r{j}{i,k};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

335

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

336 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

337 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

338 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

339

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

340 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

341 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

342 obj.T_l = T_l;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

343 obj.T_r = T_r;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

344 obj.tau_l = tau_l;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

345 obj.tau_r = tau_r;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

346

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

347 % Kroneckered norms and coefficients

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

348 I_dim = speye(dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

349 obj.RHOi_kron = kron(obj.RHOi, I_dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

350 obj.Ji_kron = kron(obj.Ji, I_dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

351 obj.Hi_kron = kron(obj.Hi, I_dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

352 obj.H_kron = kron(obj.H, I_dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

353 obj.J_kron = kron(obj.J, I_dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

354

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

355 % Misc.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

356 obj.h = g.scaling();

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

357 obj.order = order;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

358 obj.grid = g;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

359 obj.dim = dim;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

360

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

361 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

362

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

363

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

364 % Closure functions return the operators applied to the own domain to close the boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

365 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

366 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

367 % type is a cell array of strings specifying the type of boundary condition for each component.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

368 % data is a function returning the data that should be applied at the boundary.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

369 % neighbour_scheme is an instance of Scheme that should be interfaced to.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

370 % neighbour_boundary is a string specifying which boundary to interface to.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

371 function [closure, penalty] = boundary_condition(obj, boundary, type, tuning)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

372 default_arg('type',{'free','free'});

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

373 default_arg('tuning', 1.2);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

374

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

375 if ~iscell(type)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

376 type = {type, type};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

377 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

378

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

379 % j is the coordinate direction of the boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

380 j = obj.get_boundary_number(boundary);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

381 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

382

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

383 E = obj.E;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

384 Hi = obj.Hi;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

385 LAMBDA = obj.LAMBDA;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

386 MU = obj.MU;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

387 RHOi = obj.RHOi;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

388 Ji = obj.Ji;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

389

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

390 dim = obj.dim;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

391 m_tot = obj.grid.N();

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

392

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

393 % Preallocate

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

394 closure = sparse(dim*m_tot, dim*m_tot);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

395 penalty = cell(dim,1);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

396 for k = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

397 penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j));

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

398 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

399

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

400 % Loop over components that we (potentially) have different BC on

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

401 for k = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

402 switch type{k}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

403

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

404 % Dirichlet boundary condition

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

405 case {'D','d','dirichlet','Dirichlet'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

406

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

407 phi = obj.phi{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

408 h = obj.h(j);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

409 h11 = obj.H11{j}*h;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

410 gamma = obj.gamma{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

411

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

412 a_lambda = dim/h11 + 1/(h11*phi);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

413 a_mu_i = 2/(gamma*h);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

414 a_mu_ij = 2/h11 + 1/(h11*phi);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

415

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

416 d = @kroneckerDelta; % Kronecker delta

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

417 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

418 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

419 + d(i,j)* a_mu_i*MU ...

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

420 + db(i,j)*a_mu_ij*MU );

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

421

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

422 % Loop over components that Dirichlet penalties end up on

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

423 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

424 C = T{k,i};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

425 A = -d(i,k)*alpha(i,j);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

426 B = A + C;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

427 closure = closure + E{i}*RHOi*Hi*Ji*B'*e*H_gamma*(e'*E{k}' );

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

428 penalty{k} = penalty{k} - E{i}*RHOi*Hi*Ji*B'*e*H_gamma;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

429 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

430

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

431 % Free boundary condition

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

432 case {'F','f','Free','free','traction','Traction','t','T'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

433 closure = closure - E{k}*RHOi*Ji*Hi*e*H_gamma* (e'*tau{k} );

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

434 penalty{k} = penalty{k} + E{k}*RHOi*Ji*Hi*e*H_gamma;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

435

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

436 % Unknown boundary condition

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

437 otherwise

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

438 error('No such boundary condition: type = %s',type);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

439 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

440 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

441 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

442

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

443 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

444 % u denotes the solution in the own domain

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

445 % v denotes the solution in the neighbour domain

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

446 % Operators without subscripts are from the own domain.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

447 error('Not implemented');

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

448 tuning = 1.2;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

449

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

450 % j is the coordinate direction of the boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

451 j = obj.get_boundary_number(boundary);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

452 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

453

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

454 % Get boundary operators

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

455 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

456 [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

457

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

458 % Operators and quantities that correspond to the own domain only

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

459 Hi = obj.Hi;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

460 RHOi = obj.RHOi;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

461 dim = obj.dim;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

462

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

463 %--- Other operators ----

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

464 m_tot_u = obj.grid.N();

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

465 E = obj.E;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

466 LAMBDA_u = obj.LAMBDA;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

467 MU_u = obj.MU;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

468 lambda_u = e'*LAMBDA_u*e;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

469 mu_u = e'*MU_u*e;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

470

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

471 m_tot_v = neighbour_scheme.grid.N();

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

472 E_v = neighbour_scheme.E;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

473 LAMBDA_v = neighbour_scheme.LAMBDA;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

474 MU_v = neighbour_scheme.MU;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

475 lambda_v = e_v'*LAMBDA_v*e_v;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

476 mu_v = e_v'*MU_v*e_v;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

477 %-------------------------

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

478

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

479 % Borrowing constants

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

480 phi_u = obj.phi{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

481 h_u = obj.h(j);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

482 h11_u = obj.H11{j}*h_u;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

483 gamma_u = obj.gamma{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

484

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

485 phi_v = neighbour_scheme.phi{j_v};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

486 h_v = neighbour_scheme.h(j_v);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

487 h11_v = neighbour_scheme.H11{j_v}*h_v;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

488 gamma_v = neighbour_scheme.gamma{j_v};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

489

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

490 % E > sum_i 1/(2*alpha_ij)*(tau_i)^2

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

491 function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

492 th1 = h11/(2*dim);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

493 th2 = h11*phi/2;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

494 th3 = h*gamma;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

495 a1 = ( (th1 + th2)*th3*lambda + 4*th1*th2*mu ) / (2*th1*th2*th3);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

496 a2 = ( 16*(th1 + th2)*lambda*mu ) / (th1*th2*th3);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

497 alpha_ii = a1 + sqrt(a2 + a1^2);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

498

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

499 alpha_ij = mu*(2/h11 + 1/(phi*h11));

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

500 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

501

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

502 [alpha_ii_u, alpha_ij_u] = computeAlpha(phi_u,h_u,h11_u,gamma_u,lambda_u,mu_u);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

503 [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

504 sigma_ii = tuning*(alpha_ii_u + alpha_ii_v)/4;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

505 sigma_ij = tuning*(alpha_ij_u + alpha_ij_v)/4;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

506

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

507 d = @kroneckerDelta; % Kronecker delta

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

508 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

509 sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

510

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

511 % Preallocate

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

512 closure = sparse(dim*m_tot_u, dim*m_tot_u);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

513 penalty = sparse(dim*m_tot_u, dim*m_tot_v);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

514

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

515 % Loop over components that penalties end up on

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

516 for i = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

517 closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}';

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

518 penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}';

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

519

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

520 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

521 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

522

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

523 % Loop over components that we have interface conditions on

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

524 for k = 1:dim

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

525 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}';

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

526 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}';

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

527 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

528 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

529 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

530

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

531 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

532 function [j, nj] = get_boundary_number(obj, boundary)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

533

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

534 switch boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

535 case {'w','W','west','West', 'e', 'E', 'east', 'East'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

536 j = 1;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

537 case {'s','S','south','South', 'n', 'N', 'north', 'North'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

538 j = 2;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

539 otherwise

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

540 error('No such boundary: boundary = %s',boundary);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

541 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

542

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

543 switch boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

544 case {'w','W','west','West','s','S','south','South'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

545 nj = -1;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

546 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

547 nj = 1;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

548 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

549 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

550

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

551 % Returns the boundary operator op for the boundary specified by the string boundary.

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

552 % op: may be a cell array of strings

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

553 function [varargout] = get_boundary_operator(obj, op, boundary)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

554

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

555 switch boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

556 case {'w','W','west','West', 'e', 'E', 'east', 'East'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

557 j = 1;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

558 case {'s','S','south','South', 'n', 'N', 'north', 'North'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

559 j = 2;

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

560 otherwise

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

561 error('No such boundary: boundary = %s',boundary);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

562 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

563

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

564 if ~iscell(op)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

565 op = {op};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

566 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

567

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

568 for i = 1:length(op)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

569 switch op{i}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

570 case 'e'

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

571 switch boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

572 case {'w','W','west','West','s','S','south','South'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

573 varargout{i} = obj.e_l{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

574 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

575 varargout{i} = obj.e_r{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

576 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

577 case 'd'

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

578 switch boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

579 case {'w','W','west','West','s','S','south','South'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

580 varargout{i} = obj.d1_l{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

581 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

582 varargout{i} = obj.d1_r{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

583 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

584 case 'H'

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

585 switch boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

586 case {'w','W','west','West','s','S','south','South'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

587 varargout{i} = obj.H_boundary_l{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

588 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

589 varargout{i} = obj.H_boundary_r{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

590 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

591 case 'T'

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

592 switch boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

593 case {'w','W','west','West','s','S','south','South'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

594 varargout{i} = obj.T_l{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

595 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

596 varargout{i} = obj.T_r{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

597 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

598 case 'tau'

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

599 switch boundary

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

600 case {'w','W','west','West','s','S','south','South'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

601 varargout{i} = obj.tau_l{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

602 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

603 varargout{i} = obj.tau_r{j};

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

604 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

605 otherwise

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

606 error(['No such operator: operator = ' op{i}]);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

607 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

608 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

609

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

610 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

611

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

612 function N = size(obj)

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

613 N = obj.dim*prod(obj.m);

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

614 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

615 end

8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working.

Martin Almquist <malmquist@stanford.edu>

parents:

diff changeset

616 end

Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinear.m @ 739:8efc04e97da4 feature/poroelastic