sbplib: +scheme/Elastic2dVariable.m annotate

annotate +scheme/Elastic2dVariable.m @ 882:14fee299ada2 feature/poroelastic

In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.

author	Martin Almquist <malmquist@stanford.edu>
date	Fri, 26 Oct 2018 16:35:23 -0700
parents	607c631f175e
children	e30aaa4a3e09 386ef449df51

rev	line source
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	1 classdef Elastic2dVariable < scheme.Scheme
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	2
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	3 % Discretizes the elastic wave equation:
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	4 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	5 % opSet should be cell array of opSets, one per dimension. This
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	6 % is useful if we have periodic BC in one direction.
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	7
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	8 properties
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	9 m % Number of points in each direction, possibly a vector
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	10 h % Grid spacing
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	11
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	12 grid
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	13 dim
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	14
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	15 order % Order of accuracy for the approximation
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	16
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	17 % Diagonal matrices for varible coefficients
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	18 LAMBDA % Variable coefficient, related to dilation
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	19 MU % Shear modulus, variable coefficient
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	20 RHO, RHOi % Density, variable
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	21
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	22 D % Total operator
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	23 D1 % First derivatives
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	24
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	25 % Second derivatives
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	26 D2_lambda
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	27 D2_mu
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	28
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	29 % Traction operators used for BC
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	30 T_l, T_r
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	31 tau_l, tau_r
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	32
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	33 H, Hi, H_1D % Inner products
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	34 e_l, e_r
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	35 d1_l, d1_r % Normal derivatives at the boundary
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	36 E % E{i}^T picks out component i
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	37
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	38 H_boundary % Boundary inner products
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	39
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	40 % Kroneckered norms and coefficients
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	41 RHOi_kron
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	42 Hi_kron
861 607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	43
882 14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	44 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	45 theta_R % Borrowing (d1- D1)^2 from R
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	46 theta_H % First entry in norm matrix
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	47 theta_M % Borrowing d1^2 from M.
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	48
861 607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	49 % Structures used for adjoint optimization
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	50 B
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	51 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	52
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	53 methods
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	54
861 607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	55 % The coefficients can either be function handles or grid functions
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	56 function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet)
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	57 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
861 607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	58 default_arg('lambda', @(x,y) 0*x+1);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	59 default_arg('mu', @(x,y) 0*x+1);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	60 default_arg('rho', @(x,y) 0*x+1);
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	61 dim = 2;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	62
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	63 assert(isa(g, 'grid.Cartesian'))
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	64
861 607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	65 if isa(lambda, 'function_handle')
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	66 lambda = grid.evalOn(g, lambda);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	67 end
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	68 if isa(mu, 'function_handle')
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	69 mu = grid.evalOn(g, mu);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	70 end
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	71 if isa(rho, 'function_handle')
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	72 rho = grid.evalOn(g, rho);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	73 end
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	74
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	75 m = g.size();
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	76 m_tot = g.N();
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	77
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	78 h = g.scaling();
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	79 lim = g.lim;
729 aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	80 if isempty(lim)
aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	81 x = g.x;
aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	82 lim = cell(length(x),1);
aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	83 for i = 1:length(x)
aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	84 lim{i} = {min(x{i}), max(x{i})};
aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	85 end
aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	86 end
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	87
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	88 % 1D operators
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	89 ops = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	90 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	91 ops{i} = opSet{i}(m(i), lim{i}, order);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	92 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	93
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	94 % Borrowing constants
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	95 for i = 1:dim
882 14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	96 obj.theta_R{i} = h(i)*ops{i}.borrowing.R.delta_D;
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	97 obj.theta_H{i} = h(i)*ops{i}.borrowing.H11;
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	98 obj.theta_M{i} = h(i)*ops{i}.borrowing.M.d1;
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	99 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	100
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	101 I = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	102 D1 = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	103 D2 = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	104 H = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	105 Hi = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	106 e_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	107 e_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	108 d1_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	109 d1_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	110
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	111 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	112 I{i} = speye(m(i));
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	113 D1{i} = ops{i}.D1;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	114 D2{i} = ops{i}.D2;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	115 H{i} = ops{i}.H;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	116 Hi{i} = ops{i}.HI;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	117 e_l{i} = ops{i}.e_l;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	118 e_r{i} = ops{i}.e_r;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	119 d1_l{i} = ops{i}.d1_l;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	120 d1_r{i} = ops{i}.d1_r;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	121 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	122
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	123 %====== Assemble full operators ========
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	124 LAMBDA = spdiag(lambda);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	125 obj.LAMBDA = LAMBDA;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	126 MU = spdiag(mu);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	127 obj.MU = MU;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	128 RHO = spdiag(rho);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	129 obj.RHO = RHO;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	130 obj.RHOi = inv(RHO);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	131
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	132 obj.D1 = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	133 obj.D2_lambda = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	134 obj.D2_mu = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	135 obj.e_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	136 obj.e_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	137 obj.d1_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	138 obj.d1_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	139
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	140 % D1
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	141 obj.D1{1} = kron(D1{1},I{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	142 obj.D1{2} = kron(I{1},D1{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	143
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	144 % Boundary operators
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	145 obj.e_l{1} = kron(e_l{1},I{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	146 obj.e_l{2} = kron(I{1},e_l{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	147 obj.e_r{1} = kron(e_r{1},I{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	148 obj.e_r{2} = kron(I{1},e_r{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	149
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	150 obj.d1_l{1} = kron(d1_l{1},I{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	151 obj.d1_l{2} = kron(I{1},d1_l{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	152 obj.d1_r{1} = kron(d1_r{1},I{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	153 obj.d1_r{2} = kron(I{1},d1_r{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	154
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	155 % D2
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	156 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	157 obj.D2_lambda{i} = sparse(m_tot);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	158 obj.D2_mu{i} = sparse(m_tot);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	159 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	160 ind = grid.funcToMatrix(g, 1:m_tot);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	161
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	162 for i = 1:m(2)
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	163 D_lambda = D2{1}(lambda(ind(:,i)));
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	164 D_mu = D2{1}(mu(ind(:,i)));
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	165
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	166 p = ind(:,i);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	167 obj.D2_lambda{1}(p,p) = D_lambda;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	168 obj.D2_mu{1}(p,p) = D_mu;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	169 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	170
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	171 for i = 1:m(1)
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	172 D_lambda = D2{2}(lambda(ind(i,:)));
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	173 D_mu = D2{2}(mu(ind(i,:)));
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	174
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	175 p = ind(i,:);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	176 obj.D2_lambda{2}(p,p) = D_lambda;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	177 obj.D2_mu{2}(p,p) = D_mu;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	178 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	179
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	180 % Quadratures
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	181 obj.H = kron(H{1},H{2});
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	182 obj.Hi = inv(obj.H);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	183 obj.H_boundary = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	184 obj.H_boundary{1} = H{2};
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	185 obj.H_boundary{2} = H{1};
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	186 obj.H_1D = {H{1}, H{2}};
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	187
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	188 % E{i}^T picks out component i.
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	189 E = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	190 I = speye(m_tot,m_tot);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	191 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	192 e = sparse(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	193 e(i) = 1;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	194 E{i} = kron(I,e);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	195 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	196 obj.E = E;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	197
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	198 % Differentiation matrix D (without SAT)
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	199 D2_lambda = obj.D2_lambda;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	200 D2_mu = obj.D2_mu;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	201 D1 = obj.D1;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	202 D = sparse(dimm_tot,dimm_tot);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	203 d = @kroneckerDelta; % Kronecker delta
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	204 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	205 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	206 for j = 1:dim
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	207 D = D + E{i}inv(RHO)( d(i,j)D2_lambda{i}E{j}' +...
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	208 db(i,j)D1{i}LAMBDAD1{j}E{j}' ...
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	209 );
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	210 D = D + E{i}inv(RHO)( d(i,j)D2_mu{i}E{j}' +...
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	211 db(i,j)D1{j}MUD1{i}E{j}' + ...
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	212 D2_mu{j}*E{i}' ...
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	213 );
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	214 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	215 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	216 obj.D = D;
861 607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	217 %=========================================%'
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	218
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	219 % Numerical traction operators for BC.
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	220 % Because d1 =/= e0^T*D1, the numerical tractions are different
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	221 % at every boundary.
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	222 T_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	223 T_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	224 tau_l = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	225 tau_r = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	226 % tau^{j}_i = sum_k T^{j}_{ik} u_k
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	227
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	228 d1_l = obj.d1_l;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	229 d1_r = obj.d1_r;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	230 e_l = obj.e_l;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	231 e_r = obj.e_r;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	232 D1 = obj.D1;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	233
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	234 % Loop over boundaries
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	235 for j = 1:dim
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	236 T_l{j} = cell(dim,dim);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	237 T_r{j} = cell(dim,dim);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	238 tau_l{j} = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	239 tau_r{j} = cell(dim,1);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	240
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	241 % Loop over components
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	242 for i = 1:dim
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	243 tau_l{j}{i} = sparse(m_tot,dim*m_tot);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	244 tau_r{j}{i} = sparse(m_tot,dim*m_tot);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	245 for k = 1:dim
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	246 T_l{j}{i,k} = ...
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	247 -d(i,j)LAMBDA(d(i,k)e_l{k}d1_l{k}' + db(i,k)*D1{k})...
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	248 -d(j,k)MU(d(i,j)e_l{i}d1_l{i}' + db(i,j)*D1{i})...
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	249 -d(i,k)MUe_l{j}*d1_l{j}';
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	250
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	251 T_r{j}{i,k} = ...
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	252 d(i,j)LAMBDA(d(i,k)e_r{k}d1_r{k}' + db(i,k)*D1{k})...
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	253 +d(j,k)MU(d(i,j)e_r{i}d1_r{i}' + db(i,j)*D1{i})...
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	254 +d(i,k)MUe_r{j}*d1_r{j}';
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	255
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	256 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	257 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	258 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	259
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	260 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	261 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	262 obj.T_l = T_l;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	263 obj.T_r = T_r;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	264 obj.tau_l = tau_l;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	265 obj.tau_r = tau_r;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	266
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	267 % Kroneckered norms and coefficients
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	268 I_dim = speye(dim);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	269 obj.RHOi_kron = kron(obj.RHOi, I_dim);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	270 obj.Hi_kron = kron(obj.Hi, I_dim);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	271
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	272 % Misc.
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	273 obj.m = m;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	274 obj.h = h;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	275 obj.order = order;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	276 obj.grid = g;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	277 obj.dim = dim;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	278
861 607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	279 % Used for adjoint optimization
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	280 obj.B = cell(1,dim);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	281 for i = 1:dim
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	282 obj.B{i} = zeros(m(i),m(i),m(i));
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	283 for k = 1:m(i)
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	284 c = sparse(m(i),1);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	285 c(k) = 1;
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	286 [~, obj.B{i}(:,:,k)] = ops{i}.D2(c);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	287 end
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	288 end
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	289
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	290 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	291
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	292
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	293 % Closure functions return the operators applied to the own domain to close the boundary
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	294 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	295 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	296 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	297 % on the first component.
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	298 % data is a function returning the data that should be applied at the boundary.
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	299 % neighbour_scheme is an instance of Scheme that should be interfaced to.
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	300 % neighbour_boundary is a string specifying which boundary to interface to.
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	301 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
738 aa4ef495f1fd Set Dirichlet tuning as parameter in BC function in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 734 diff changeset	302 default_arg('tuning', 1.2);
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	303
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	304 assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	305 comp = bc{1};
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	306 type = bc{2};
729 aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	307
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	308 % j is the coordinate direction of the boundary
726 37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	309 j = obj.get_boundary_number(boundary);
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	310 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	311
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	312 E = obj.E;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	313 Hi = obj.Hi;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	314 LAMBDA = obj.LAMBDA;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	315 MU = obj.MU;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	316 RHOi = obj.RHOi;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	317
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	318 dim = obj.dim;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	319 m_tot = obj.grid.N();
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	320
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	321 % Preallocate
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	322 closure = sparse(dimm_tot, dimm_tot);
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	323 penalty = sparse(dim*m_tot, m_tot/obj.m(j));
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	324
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	325 k = comp;
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	326 switch type
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	327
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	328 % Dirichlet boundary condition
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	329 case {'D','d','dirichlet','Dirichlet'}
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	330
882 14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	331 theta_R = obj.theta_R{j};
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	332 theta_H = obj.theta_H{j};
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	333 theta_M = obj.theta_M{j};
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	334
882 14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	335 a_lambda = dim/theta_H + 1/theta_R;
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	336 a_mu_i = 2/theta_M;
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	337 a_mu_ij = 2/theta_H + 1/theta_R;
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	338
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	339 d = @kroneckerDelta; % Kronecker delta
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	340 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	341 alpha = @(i,j) tuning( d(i,j) a_lambda*LAMBDA ...
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	342 + d(i,j)* a_mu_i*MU ...
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	343 + db(i,j)a_mu_ijMU );
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	344
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	345 % Loop over components that Dirichlet penalties end up on
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	346 for i = 1:dim
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	347 C = T{k,i};
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	348 A = -d(i,k)*alpha(i,j);
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	349 B = A + C;
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	350 closure = closure + E{i}RHOiHiB'eH_gamma(e'*E{k}' );
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	351 penalty = penalty - E{i}RHOiHiB'e*H_gamma;
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	352 end
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	353
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	354 % Free boundary condition
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	355 case {'F','f','Free','free','traction','Traction','t','T'}
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	356 closure = closure - E{k}RHOiHieH_gamma* (e'*tau{k} );
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	357 penalty = penalty + E{k}RHOiHieH_gamma;
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	358
795 1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	359 % Unknown boundary condition
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	360 otherwise
1f6b2fb69225 Revert bcSetup and update bc functions in elastic schemes to be compatible. Martin Almquist <malmquist@stanford.edu> parents: 738 diff changeset	361 error('No such boundary condition: type = %s',type);
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	362 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	363 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	364
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	365 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	366 % u denotes the solution in the own domain
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	367 % v denotes the solution in the neighbour domain
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	368 % Operators without subscripts are from the own domain.
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	369 tuning = 1.2;
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	370
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	371 % j is the coordinate direction of the boundary
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	372 j = obj.get_boundary_number(boundary);
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	373 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	374
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	375 % Get boundary operators
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	376 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	377 [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	378
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	379 % Operators and quantities that correspond to the own domain only
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	380 Hi = obj.Hi;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	381 RHOi = obj.RHOi;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	382 dim = obj.dim;
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	383
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	384 %--- Other operators ----
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	385 m_tot_u = obj.grid.N();
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	386 E = obj.E;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	387 LAMBDA_u = obj.LAMBDA;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	388 MU_u = obj.MU;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	389 lambda_u = e'LAMBDA_ue;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	390 mu_u = e'MU_ue;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	391
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	392 m_tot_v = neighbour_scheme.grid.N();
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	393 E_v = neighbour_scheme.E;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	394 LAMBDA_v = neighbour_scheme.LAMBDA;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	395 MU_v = neighbour_scheme.MU;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	396 lambda_v = e_v'LAMBDA_ve_v;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	397 mu_v = e_v'MU_ve_v;
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	398 %-------------------------
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	399
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	400 % Borrowing constants
882 14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	401 theta_R_u = obj.theta_R{j};
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	402 theta_H_u = obj.theta_H{j};
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	403 theta_M_u = obj.theta_M{j};
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	404
882 14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	405 theta_R_v = neighbour_scheme.theta_R{j_v};
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	406 theta_H_v = neighbour_scheme.theta_H{j_v};
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	407 theta_M_v = neighbour_scheme.theta_M{j_v};
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	408
882 14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	409 function [alpha_ii, alpha_ij] = computeAlpha(th_R, th_H, th_M, lambda, mu)
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	410 alpha_ii = dimlambda/(4th_H) + lambda/(4th_R) + mu/(2th_M);
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	411 alpha_ij = mu/(2th_H) + mu/(4th_R);
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	412 end
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	413
882 14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	414 [alpha_ii_u, alpha_ij_u] = computeAlpha(theta_R_u, theta_H_u, theta_M_u, lambda_u, mu_u);
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	415 [alpha_ii_v, alpha_ij_v] = computeAlpha(theta_R_v, theta_H_v, theta_M_v, lambda_v, mu_v);
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	416 sigma_ii = tuning*(alpha_ii_u + alpha_ii_v);
14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	417 sigma_ij = tuning*(alpha_ij_u + alpha_ij_v);
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	418
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	419 d = @kroneckerDelta; % Kronecker delta
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	420 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	421 sigma = @(i,j) tuning(d(i,j)sigma_ii + db(i,j)*sigma_ij);
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	422
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	423 % Preallocate
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	424 closure = sparse(dimm_tot_u, dimm_tot_u);
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	425 penalty = sparse(dimm_tot_u, dimm_tot_v);
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	426
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	427 % Loop over components that penalties end up on
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	428 for i = 1:dim
729 aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	429 closure = closure - E{i}RHOiHiesigma(i,j)H_gammae'*E{i}';
aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	430 penalty = penalty + E{i}RHOiHiesigma(i,j)H_gammae_v'*E_v{i}';
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	431
729 aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	432 closure = closure - 1/2E{i}RHOiHieH_gammae'*tau{i};
aa8cf3851de8 Update multiblock.DiffOp to work for systems. Martin Almquist <malmquist@stanford.edu> parents: 727 diff changeset	433 penalty = penalty - 1/2E{i}RHOiHieH_gammae_v'*tau_v{i};
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	434
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	435 % Loop over components that we have interface conditions on
6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	436 for k = 1:dim
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	437 closure = closure + 1/2E{i}RHOiHiT{k,i}'eH_gammae'E{k}';
5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	438 penalty = penalty - 1/2E{i}RHOiHiT{k,i}'eH_gammae_v'E_v{k}';
5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	439 end
727 6d5953fc090e First implementation of elastic interface Martin Almquist <malmquist@stanford.edu> parents: 726 diff changeset	440 end
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	441 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	442
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	443 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	444 function [j, nj] = get_boundary_number(obj, boundary)
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	445
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	446 switch boundary
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	447 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	448 j = 1;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	449 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	450 j = 2;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	451 otherwise
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	452 error('No such boundary: boundary = %s',boundary);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	453 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	454
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	455 switch boundary
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	456 case {'w','W','west','West','s','S','south','South'}
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	457 nj = -1;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	458 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	459 nj = 1;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	460 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	461 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	462
726 37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	463 % Returns the boundary operator op for the boundary specified by the string boundary.
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	464 % op: may be a cell array of strings
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	465 function [varargout] = get_boundary_operator(obj, op, boundary)
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	466
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	467 switch boundary
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	468 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	469 j = 1;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	470 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	471 j = 2;
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	472 otherwise
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	473 error('No such boundary: boundary = %s',boundary);
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	474 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	475
726 37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	476 if ~iscell(op)
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	477 op = {op};
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	478 end
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	479
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	480 for i = 1:length(op)
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	481 switch op{i}
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	482 case 'e'
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	483 switch boundary
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	484 case {'w','W','west','West','s','S','south','South'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	485 varargout{i} = obj.e_l{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	486 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	487 varargout{i} = obj.e_r{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	488 end
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	489 case 'd'
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	490 switch boundary
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	491 case {'w','W','west','West','s','S','south','South'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	492 varargout{i} = obj.d1_l{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	493 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	494 varargout{i} = obj.d1_r{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	495 end
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	496 case 'H'
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	497 varargout{i} = obj.H_boundary{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	498 case 'T'
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	499 switch boundary
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	500 case {'w','W','west','West','s','S','south','South'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	501 varargout{i} = obj.T_l{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	502 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	503 varargout{i} = obj.T_r{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	504 end
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	505 case 'tau'
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	506 switch boundary
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	507 case {'w','W','west','West','s','S','south','South'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	508 varargout{i} = obj.tau_l{j};
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	509 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	510 varargout{i} = obj.tau_r{j};
855 5751262b323b Add 1D quadrature matrices as property in Elastic2dVariable. Martin Almquist <malmquist@stanford.edu> parents: 795 diff changeset	511 end
861 607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	512 case 'alpha'
882 14fee299ada2 In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok. Martin Almquist <malmquist@stanford.edu> parents: 861 diff changeset	513 % alpha = alpha(i,j) is the penalty strength for displacement BC.
861 607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	514 tuning = 1.2;
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	515 LAMBDA = obj.LAMBDA;
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	516 MU = obj.MU;
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	517
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	518 phi = obj.phi{j};
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	519 h = obj.h(j);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	520 h11 = obj.H11{j}*h;
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	521 gamma = obj.gamma{j};
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	522 dim = obj.dim;
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	523
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	524 a_lambda = dim/h11 + 1/(h11*phi);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	525 a_mu_i = 2/(gamma*h);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	526 a_mu_ij = 2/h11 + 1/(h11*phi);
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	527
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	528 d = @kroneckerDelta; % Kronecker delta
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	529 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	530 alpha = @(i,k) d(i,k)tuning( d(i,j)* a_lambda*LAMBDA ...
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	531 + d(i,j)* a_mu_i*MU ...
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	532 + db(i,j)a_mu_ijMU );
607c631f175e Small changes to Elastic2dVariable to facilitate adjoing gradient computation. Martin Almquist <malmquist@stanford.edu> parents: 855 diff changeset	533 varargout{i} = alpha;
726 37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	534 otherwise
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	535 error(['No such operator: operator = ' op{i}]);
37d5d69b1a0d Refactor BC code in Elastic2dVariable Martin Almquist <malmquist@stanford.edu> parents: 689 diff changeset	536 end
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	537 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	538
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	539 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	540
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	541 function N = size(obj)
734 eebe24a636c7 Make Elastic2dVariable.size account for components. Martin Almquist <malmquist@stanford.edu> parents: 730 diff changeset	542 N = obj.dim*prod(obj.m);
687 e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	543 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	544 end
e8fc3aa1faf6 Rename elastic scheme. Martin Almquist <malmquist@stanford.edu> parents: diff changeset	545 end

Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dVariable.m @ 882:14fee299ada2 feature/poroelastic