sbplib: +scheme/Elastic2dCurvilinear.m annotate

annotate +scheme/Elastic2dCurvilinear.m @ 740:f4e2a6a2df08 feature/poroelastic

Make computation of metric derivatives work for periodic.

author	Martin Almquist <malmquist@stanford.edu>
date	Mon, 07 May 2018 15:55:35 -0700
parents	8efc04e97da4
children	08f3ffe63f48

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
740 f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	141 % Use non-periodic difference operators for metric even if opSet is periodic.
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	142 xmax = max(ops{1}.x);
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	143 ymax = max(ops{2}.x);
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	144 opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order);
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	145 opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order);
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	146 D1Metric{1} = kron(opSetMetric{1}.D1, I{2});
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	147 D1Metric{2} = kron(I{1}, opSetMetric{2}.D1);
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	148
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	149 x_xi = D1Metric{1}*x;
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	150 x_eta = D1Metric{2}*x;
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	151 y_xi = D1Metric{1}*y;
f4e2a6a2df08 Make computation of metric derivatives work for periodic. Martin Almquist <malmquist@stanford.edu> parents: 739 diff changeset	152 y_eta = D1Metric{2}*y;
739 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 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	155
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	156 b = cell(dim,dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	157 b{1,1} = y_eta./J;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	158 b{1,2} = -x_eta./J;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	159 b{2,1} = -y_xi./J;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	160 b{2,2} = x_xi./J;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	161
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	162 % Scale factors for boundary integrals
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	163 beta = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	164 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	165 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	166
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	167 J = spdiag(J);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	168 Ji = inv(J);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	169 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	170 beta{i} = spdiag(beta{i});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	171 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	172 b{i,j} = spdiag(b{i,j});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	173 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	174 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	175 obj.J = J;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	176 obj.Ji = Ji;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	177 obj.b = b;
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
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	180 % Boundary operators
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	181 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	182 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	183 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	184 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	185
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	186 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	187 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	188 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	189 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	190
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	191 % D2
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	192 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	193 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	194 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	195 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	196 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	197 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	198 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	199 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	200 ind = grid.funcToMatrix(g, 1:m_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	201
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	202 % x-dir
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	203 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	204 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	205 for k = 1
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 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	208 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	209
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	210 for col = 1:m(2)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	211 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	212 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	213
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	214 p = ind(:,col);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	215 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	216 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	217 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	218
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	219 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	220 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	221 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	222
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	223 % y-dir
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	224 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	225 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	226 for k = 2
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 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	229 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	230
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	231 for row = 1:m(1)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	232 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	233 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	234
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	235 p = ind(row,:);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	236 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	237 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	238 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	239
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	240 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	241 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	242 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	243
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	244 % Quadratures
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	245 obj.H = kron(H{1},H{2});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	246 obj.Hi = inv(obj.H);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	247 obj.H_boundary_l = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	248 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	249 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	250 obj.H_boundary_r = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	251 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	252 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	253
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	254 % E{i}^T picks out component i.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	255 E = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	256 I = speye(m_tot,m_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	257 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	258 e = sparse(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	259 e(i) = 1;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	260 E{i} = kron(I,e);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	261 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	262 obj.E = E;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	263
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	264 % Differentiation matrix D (without SAT)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	265 D2_lambda = obj.D2_lambda;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	266 D2_mu = obj.D2_mu;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	267 D1 = obj.D1;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	268 D = sparse(dimm_tot,dimm_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	269 d = @kroneckerDelta; % Kronecker delta
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	270 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	271 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	272 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	273 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	274 for l = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	275 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	276 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	277 );
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 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	280 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	281 );
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	282
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	283 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	284 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	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 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	288 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	289 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	290 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	291 obj.D = D;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	292 %=========================================%
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	293
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	294 % Numerical traction operators for BC.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	295 % 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	296 % at every boundary.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	297 T_l = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	298 T_r = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	299 tau_l = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	300 tau_r = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	301 % 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	302
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	303 d1_l = obj.d1_l;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	304 d1_r = obj.d1_r;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	305 e_l = obj.e_l;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	306 e_r = obj.e_r;
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 boundaries
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	309 for j = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	310 T_l{j} = cell(dim,dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	311 T_r{j} = cell(dim,dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	312 tau_l{j} = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	313 tau_r{j} = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	314
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	315 % Loop over components
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	316 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	317 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	318 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	319
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	320 % 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	321 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	322
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	323 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	324 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	325
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	326 for m = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	327 for l = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	328 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	329 -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	330 -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	331 -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	332
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	333 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	334 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	335 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	336 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	337 end
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 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	341 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	342
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	343 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	344 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	345 end
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 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	348 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	349 obj.T_l = T_l;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	350 obj.T_r = T_r;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	351 obj.tau_l = tau_l;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	352 obj.tau_r = tau_r;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	353
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	354 % Kroneckered norms and coefficients
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	355 I_dim = speye(dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	356 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	357 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	358 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	359 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	360 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	361
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	362 % Misc.
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	363 obj.h = g.scaling();
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	364 obj.order = order;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	365 obj.grid = g;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	366 obj.dim = dim;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	367
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	368 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	369
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	370
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	371 % 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	372 % 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	373 % 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	374 % 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	375 % 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	376 % 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	377 % 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	378 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	379 default_arg('type',{'free','free'});
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	380 default_arg('tuning', 1.2);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	381
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	382 if ~iscell(type)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	383 type = {type, type};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	384 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	385
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	386 % 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	387 j = obj.get_boundary_number(boundary);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	388 [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	389
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	390 E = obj.E;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	391 Hi = obj.Hi;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	392 LAMBDA = obj.LAMBDA;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	393 MU = obj.MU;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	394 RHOi = obj.RHOi;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	395 Ji = obj.Ji;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	396
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	397 dim = obj.dim;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	398 m_tot = obj.grid.N();
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 % Preallocate
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	401 closure = sparse(dimm_tot, dimm_tot);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	402 penalty = cell(dim,1);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	403 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	404 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	405 end
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 % 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	408 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	409 switch type{k}
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	410
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	411 % Dirichlet boundary condition
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	412 case {'D','d','dirichlet','Dirichlet'}
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	413
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	414 phi = obj.phi{j};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	415 h = obj.h(j);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	416 h11 = obj.H11{j}*h;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	417 gamma = obj.gamma{j};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	418
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	419 a_lambda = dim/h11 + 1/(h11*phi);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	420 a_mu_i = 2/(gamma*h);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	421 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	422
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	423 d = @kroneckerDelta; % Kronecker delta
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	424 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	425 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	426 + d(i,j)* a_mu_i*MU ...
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	427 + db(i,j)a_mu_ijMU );
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	428
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	429 % 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	430 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	431 C = T{k,i};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	432 A = -d(i,k)*alpha(i,j);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	433 B = A + C;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	434 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	435 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	436 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	437
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	438 % Free boundary condition
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	439 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	440 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	441 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	442
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	443 % Unknown boundary condition
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	444 otherwise
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	445 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	446 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	447 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	448 end
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 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	451 % 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	452 % 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	453 % 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	454 error('Not implemented');
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	455 tuning = 1.2;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	456
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	457 % 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	458 j = obj.get_boundary_number(boundary);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	459 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	460
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	461 % Get boundary operators
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	462 [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	463 [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	464
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	465 % 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	466 Hi = obj.Hi;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	467 RHOi = obj.RHOi;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	468 dim = obj.dim;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	469
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	470 %--- Other operators ----
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	471 m_tot_u = obj.grid.N();
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	472 E = obj.E;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	473 LAMBDA_u = obj.LAMBDA;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	474 MU_u = obj.MU;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	475 lambda_u = e'LAMBDA_ue;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	476 mu_u = e'MU_ue;
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 m_tot_v = neighbour_scheme.grid.N();
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	479 E_v = neighbour_scheme.E;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	480 LAMBDA_v = neighbour_scheme.LAMBDA;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	481 MU_v = neighbour_scheme.MU;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	482 lambda_v = e_v'LAMBDA_ve_v;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	483 mu_v = e_v'MU_ve_v;
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
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	486 % Borrowing constants
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	487 phi_u = obj.phi{j};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	488 h_u = obj.h(j);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	489 h11_u = obj.H11{j}*h_u;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	490 gamma_u = obj.gamma{j};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	491
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	492 phi_v = neighbour_scheme.phi{j_v};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	493 h_v = neighbour_scheme.h(j_v);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	494 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	495 gamma_v = neighbour_scheme.gamma{j_v};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	496
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	497 % 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	498 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	499 th1 = h11/(2*dim);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	500 th2 = h11*phi/2;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	501 th3 = h*gamma;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	502 a1 = ( (th1 + th2)th3lambda + 4th1th2mu ) / (2th1th2th3);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	503 a2 = ( 16(th1 + th2)lambdamu ) / (th1th2*th3);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	504 alpha_ii = a1 + sqrt(a2 + a1^2);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	505
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	506 alpha_ij = mu(2/h11 + 1/(phih11));
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	507 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	508
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	509 [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	510 [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	511 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	512 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	513
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	514 d = @kroneckerDelta; % Kronecker delta
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	515 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	516 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	517
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	518 % Preallocate
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	519 closure = sparse(dimm_tot_u, dimm_tot_u);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	520 penalty = sparse(dimm_tot_u, dimm_tot_v);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	521
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	522 % 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	523 for i = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	524 closure = closure - E{i}RHOiHiesigma(i,j)H_gammae'*E{i}';
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	525 penalty = penalty + E{i}RHOiHiesigma(i,j)H_gammae_v'*E_v{i}';
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	526
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	527 closure = closure - 1/2E{i}RHOiHieH_gammae'*tau{i};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	528 penalty = penalty - 1/2E{i}RHOiHieH_gammae_v'*tau_v{i};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	529
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	530 % 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	531 for k = 1:dim
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	532 closure = closure + 1/2E{i}RHOiHiT{k,i}'eH_gammae'E{k}';
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	533 penalty = penalty - 1/2E{i}RHOiHiT{k,i}'eH_gammae_v'E_v{k}';
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	534 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	535 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	536 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	537
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	538 % 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	539 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	540
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	541 switch boundary
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	542 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	543 j = 1;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	544 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	545 j = 2;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	546 otherwise
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	547 error('No such boundary: boundary = %s',boundary);
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
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	550 switch boundary
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	551 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	552 nj = -1;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	553 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	554 nj = 1;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	555 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	556 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	557
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	558 % 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	559 % 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	560 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	561
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	562 switch boundary
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	563 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	564 j = 1;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	565 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	566 j = 2;
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	567 otherwise
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	568 error('No such boundary: boundary = %s',boundary);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	569 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	570
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	571 if ~iscell(op)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	572 op = {op};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	573 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	574
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	575 for i = 1:length(op)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	576 switch op{i}
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	577 case 'e'
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.e_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.e_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 'd'
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.d1_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.d1_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 'H'
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.H_boundary_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.H_boundary_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 'T'
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.T_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.T_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 case 'tau'
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	606 switch boundary
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	607 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	608 varargout{i} = obj.tau_l{j};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	609 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	610 varargout{i} = obj.tau_r{j};
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	611 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	612 otherwise
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	613 error(['No such operator: operator = ' op{i}]);
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
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	617 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	618
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	619 function N = size(obj)
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	620 N = obj.dim*prod(obj.m);
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	621 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	622 end
8efc04e97da4 Add Elastic curvilinear. Traction and Dirichlet BC working. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	623 end

Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dCurvilinear.m @ 740:f4e2a6a2df08 feature/poroelastic