sbplib: +scheme/Elastic2dVariableAnisotropic.m annotate

annotate +scheme/Elastic2dVariableAnisotropic.m @ 1253:89dad61cad22 feature/poroelastic

Make Elastic2dVariable faster and more memory efficient

author	Martin Almquist <malmquist@stanford.edu>
date	Tue, 04 Feb 2020 10:15:42 -0800
parents	8fc2f9a4c882
children	7fbee3de0ab2

rev	line source
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	1 classdef Elastic2dVariableAnisotropic < scheme.Scheme
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	2
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	3 % Discretizes the elastic wave equation:
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	5 % opSet should be cell array of opSets, one per dimension. This
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	6 % is useful if we have periodic BC in one direction.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	7 % Assumes fully compatible operators
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	8
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	9 properties
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	10 m % Number of points in each direction, possibly a vector
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	11 h % Grid spacing
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	12
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	13 grid
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	14 dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	15
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	16 order % Order of accuracy for the approximation
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	17
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	18 % Diagonal matrices for variable coefficients
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	19 RHO, RHOi, RHOi_kron % Density
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	20 C % Elastic stiffness tensor
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	21
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	22 D % Total operator
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	23 D1 % First derivatives
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	24 % D2 % Second derivatives
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	25
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	26 % Boundary operators in cell format, used for BC
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	27 T_w, T_e, T_s, T_n
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	28
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	29 % Traction operators
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	30 tau_w, tau_e, tau_s, tau_n % Return vector field
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	31 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	32 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	33
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	34 % Inner products
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	35 H, Hi, Hi_kron, H_1D
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	36
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	37 % Boundary inner products (for scalar field)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	38 H_w, H_e, H_s, H_n
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	39
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	40 % Boundary restriction operators
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	41 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	42 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	43 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	44 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	45
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	46 % E{i}^T picks out component i
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	47 E
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	48
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	49 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	50 h11 % First entry in norm matrix
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	51
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	52 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	53
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	54 methods
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	55
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	56 % The coefficients can either be function handles or grid functions
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	57 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	58 function obj = Elastic2dVariableAnisotropic(g, order, rho, C, opSet, optFlag)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	59 default_arg('rho', @(x,y) 0*x+1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	60 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	61 default_arg('optFlag', false);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	62 dim = 2;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	63
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	64 C_default = cell(dim,dim,dim,dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	65 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	66 for j = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	67 for k = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	68 for l = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	69 C_default{i,j,k,l} = @(x,y) 0*x + 1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	70 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	71 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	72 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	73 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	74 default_arg('C', C_default);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	75 assert(isa(g, 'grid.Cartesian'))
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	76
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	77 if isa(rho, 'function_handle')
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	78 rho = grid.evalOn(g, rho);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	79 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	80
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	81 C_mat = cell(dim,dim,dim,dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	82 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	83 for j = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	84 for k = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	85 for l = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	86 if isa(C{i,j,k,l}, 'function_handle')
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	87 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	88 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	89 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	90 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	91 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	92 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	93 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	94 obj.C = C_mat;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	95
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	96 m = g.size();
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	97 m_tot = g.N();
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	98 lim = g.lim;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	99 if isempty(lim)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	100 x = g.x;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	101 lim = cell(length(x),1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	102 for i = 1:length(x)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	103 lim{i} = {min(x{i}), max(x{i})};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	104 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	105 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	106
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	107 % 1D operators
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	108 ops = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	109 h = zeros(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	110 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	111 ops{i} = opSet{i}(m(i), lim{i}, order);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	112 h(i) = ops{i}.h;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	113 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	114
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	115 % Borrowing constants
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	116 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	117 obj.h11{i} = h(i)*ops{i}.borrowing.H11;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	118 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	119
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	120 I = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	121 D1 = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	122 D2 = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	123 H = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	124 Hi = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	125 e_0 = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	126 e_m = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	127 d1_0 = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	128 d1_m = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	129
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	130 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	131 I{i} = speye(m(i));
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	132 D1{i} = ops{i}.D1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	133 D2{i} = ops{i}.D2;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	134 H{i} = ops{i}.H;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	135 Hi{i} = ops{i}.HI;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	136 e_0{i} = ops{i}.e_l;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	137 e_m{i} = ops{i}.e_r;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	138 d1_0{i} = ops{i}.d1_l;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	139 d1_m{i} = ops{i}.d1_r;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	140 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	141
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	142 %====== Assemble full operators ========
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	143 I_dim = speye(dim, dim);
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	144 RHO = spdiag(rho);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	145 obj.RHO = RHO;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	146 obj.RHOi = inv(RHO);
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	147 obj.RHOi_kron = kron(obj.RHOi, I_dim);
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	148
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	149 obj.D1 = cell(dim,1);
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	150 D2_temp = cell(dim,dim,dim);
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	151
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	152 % D1
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	153 obj.D1{1} = kron(D1{1},I{2});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	154 obj.D1{2} = kron(I{1},D1{2});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	155
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	156 % Boundary restriction operators
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	157 e_l = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	158 e_r = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	159 e_l{1} = kron(e_0{1}, I{2});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	160 e_l{2} = kron(I{1}, e_0{2});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	161 e_r{1} = kron(e_m{1}, I{2});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	162 e_r{2} = kron(I{1}, e_m{2});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	163
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	164 e_scalar_w = e_l{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	165 e_scalar_e = e_r{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	166 e_scalar_s = e_l{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	167 e_scalar_n = e_r{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	168
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	169 e_w = kron(e_scalar_w, I_dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	170 e_e = kron(e_scalar_e, I_dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	171 e_s = kron(e_scalar_s, I_dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	172 e_n = kron(e_scalar_n, I_dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	173
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	174 % E{i}^T picks out component i.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	175 E = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	176 I = speye(m_tot,m_tot);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	177 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	178 e = sparse(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	179 e(i) = 1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	180 E{i} = kron(I,e);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	181 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	182 obj.E = E;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	183
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	184 e1_w = (e_scalar_w'*E{1}')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	185 e1_e = (e_scalar_e'*E{1}')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	186 e1_s = (e_scalar_s'*E{1}')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	187 e1_n = (e_scalar_n'*E{1}')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	188
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	189 e2_w = (e_scalar_w'*E{2}')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	190 e2_e = (e_scalar_e'*E{2}')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	191 e2_s = (e_scalar_s'*E{2}')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	192 e2_n = (e_scalar_n'*E{2}')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	193
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	194
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	195 % D2
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	196 switch order
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	197 case 2
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	198 width = 3;
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	199 case 4
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	200 width = 5;
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	201 case 6
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	202 width = 7;
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	203 end
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	204 for j = 1:dim
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	205 for k = 1:dim
6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	206 for l = 1:dim
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	207 D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot);
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	208 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	209 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	210 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	211 ind = grid.funcToMatrix(g, 1:m_tot);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	212
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	213 k = 1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	214 for r = 1:m(2)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	215 p = ind(:,r);
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	216 for j = 1:dim
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	217 for l = 1:dim
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	218 coeff = C{k,j,k,l};
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	219 D_kk = D2{1}(coeff(p));
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	220 D2_temp{j,k,l}(p,p) = D_kk;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	221 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	222 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	223 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	224
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	225 k = 2;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	226 for r = 1:m(1)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	227 p = ind(r,:);
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	228 for j = 1:dim
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	229 for l = 1:dim
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	230 coeff = C{k,j,k,l};
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	231 D_kk = D2{2}(coeff(p));
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	232 D2_temp{j,k,l}(p,p) = D_kk;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	233 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	234 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	235 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	236
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	237 % Quadratures
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	238 obj.H = kron(H{1},H{2});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	239 obj.Hi = inv(obj.H);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	240 obj.H_w = H{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	241 obj.H_e = H{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	242 obj.H_s = H{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	243 obj.H_n = H{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	244 obj.H_1D = {H{1}, H{2}};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	245
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	246 % Differentiation matrix D (without SAT)
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	247 D2_temp;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	248 D1 = obj.D1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	249 D = sparse(dimm_tot,dimm_tot);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	250 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	251 for j = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	252 for k = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	253 for l = 1:dim
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	254 if i == k
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	255 D = D + E{j}D2_temp{j,k,l}E{l}';
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	256 D2_temp{j,k,l} = [];
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	257 else
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	258 D = D + E{j}D1{i}C_mat{i,j,k,l}D1{k}E{l}';
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	259 end
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	260 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	261 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	262 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	263 end
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	264 clear D2_temp;
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	265 D = obj.RHOi_kron*D;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	266 obj.D = D;
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	267 clear D;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	268 %=========================================%'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	269
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	270 % Numerical traction operators for BC.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	271 %
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	272 % Formula at boundary j: % tau^{j}_i = sum_l T^{j}_{il} u_l
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	273 %
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	274 T_l = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	275 T_r = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	276 tau_l = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	277 tau_r = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	278
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	279 D1 = obj.D1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	280
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	281 % Boundary j
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	282 for j = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	283 T_l{j} = cell(dim,dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	284 T_r{j} = cell(dim,dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	285 tau_l{j} = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	286 tau_r{j} = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	287
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	288 [~, n_l] = size(e_l{j});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	289 [~, n_r] = size(e_r{j});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	290
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	291 % Traction component i
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	292 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	293 tau_l{j}{i} = sparse(dim*m_tot, n_l);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	294 tau_r{j}{i} = sparse(dim*m_tot, n_r);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	295
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	296 % Displacement component l
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	297 for l = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	298 T_l{j}{i,l} = sparse(m_tot, n_l);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	299 T_r{j}{i,l} = sparse(m_tot, n_r);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	300
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	301 % Derivative direction k
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	302 for k = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	303 T_l{j}{i,l} = T_l{j}{i,l} ...
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	304 - (e_l{j}'C_mat{j,i,k,l}D1{k})';
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	305 T_r{j}{i,l} = T_r{j}{i,l} ...
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	306 + (e_r{j}'C_mat{j,i,k,l}D1{k})';
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	307 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	308 tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	309 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	310 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	311 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	312 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	313
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	314 % Traction tensors, T_ij
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	315 obj.T_w = T_l{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	316 obj.T_e = T_r{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	317 obj.T_s = T_l{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	318 obj.T_n = T_r{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	319
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	320 % Restriction operators
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	321 obj.e_w = e_w;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	322 obj.e_e = e_e;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	323 obj.e_s = e_s;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	324 obj.e_n = e_n;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	325
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	326 obj.e1_w = e1_w;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	327 obj.e1_e = e1_e;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	328 obj.e1_s = e1_s;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	329 obj.e1_n = e1_n;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	330
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	331 obj.e2_w = e2_w;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	332 obj.e2_e = e2_e;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	333 obj.e2_s = e2_s;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	334 obj.e2_n = e2_n;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	335
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	336 obj.e_scalar_w = e_scalar_w;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	337 obj.e_scalar_e = e_scalar_e;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	338 obj.e_scalar_s = e_scalar_s;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	339 obj.e_scalar_n = e_scalar_n;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	340
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	341 % First component of traction
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	342 obj.tau1_w = tau_l{1}{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	343 obj.tau1_e = tau_r{1}{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	344 obj.tau1_s = tau_l{2}{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	345 obj.tau1_n = tau_r{2}{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	346
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	347 % Second component of traction
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	348 obj.tau2_w = tau_l{1}{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	349 obj.tau2_e = tau_r{1}{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	350 obj.tau2_s = tau_l{2}{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	351 obj.tau2_n = tau_r{2}{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	352
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	353 % Traction vectors
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	354 obj.tau_w = (e_w'e1_wobj.tau1_w')' + (e_w'e2_wobj.tau2_w')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	355 obj.tau_e = (e_e'e1_eobj.tau1_e')' + (e_e'e2_eobj.tau2_e')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	356 obj.tau_s = (e_s'e1_sobj.tau1_s')' + (e_s'e2_sobj.tau2_s')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	357 obj.tau_n = (e_n'e1_nobj.tau1_n')' + (e_n'e2_nobj.tau2_n')';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	358
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	359 % Kroneckered norms and coefficients
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	360 obj.Hi_kron = kron(obj.Hi, I_dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	361
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	362 % Misc.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	363 obj.m = m;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	364 obj.h = h;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	365 obj.order = order;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	366 obj.grid = g;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	367 obj.dim = dim;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	368
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	369 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	370
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	371
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	372 % Closure functions return the operators applied to the own domain to close the boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	373 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	374 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	375 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	376 % on the first component. Can also be e.g.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	377 % {'normal', 'd'} or {'tangential', 't'} for conditions on
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	378 % tangential/normal component.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	379 % data is a function returning the data that should be applied at the boundary.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	380 % neighbour_scheme is an instance of Scheme that should be interfaced to.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	381 % neighbour_boundary is a string specifying which boundary to interface to.
1203 25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	382
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	383 % For displacement bc:
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	384 % bc = {comp, 'd', dComps},
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	385 % where
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	386 % dComps = vector of components with displacement BC. Default: 1:dim.
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	387 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	388 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	389 default_arg('tuning', 1.0);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	390
1202 31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	391 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	392 comp = bc{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	393 type = bc{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	394 if ischar(comp)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	395 comp = obj.getComponent(comp, boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	396 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	397
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	398 e = obj.getBoundaryOperatorForScalarField('e', boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	399 tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	400 T = obj.getBoundaryTractionOperator(boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	401 h11 = obj.getBorrowing(boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	402 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	403 nu = obj.getNormal(boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	404
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	405 E = obj.E;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	406 Hi = obj.Hi;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	407 RHOi = obj.RHOi;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	408 C = obj.C;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	409
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	410 dim = obj.dim;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	411 m_tot = obj.grid.N();
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	412
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	413 % Preallocate
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	414 [~, col] = size(tau);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	415 closure = sparse(dimm_tot, dimm_tot);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	416 penalty = sparse(dim*m_tot, col);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	417
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	418 j = comp;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	419 switch type
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	420
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	421 % Dirichlet boundary condition
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	422 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	423
1202 31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	424 if numel(bc) >= 3
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	425 dComps = bc{3};
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	426 else
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	427 dComps = 1:dim;
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	428 end
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	429
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	430 % Loops over components that Dirichlet penalties end up on
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	431 % Y: symmetrizing part of penalty
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	432 % Z: symmetric part of penalty
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	433 % X = Y + Z.
1202 31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	434
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	435 % Nonsymmetric part goes on all components to
1203 25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	436 % yield traction in discrete energy rate
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	437 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	438 Y = T{j,i}';
1202 31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	439 X = e*Y;
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	440 closure = closure + E{i}RHOiHiX'eH_gamma(e'*E{j}' );
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	441 penalty = penalty - E{i}RHOiHiX'e*H_gamma;
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	442 end
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	443
1203 25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	444 % Symmetric part only required on components with displacement BC.
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	445 % (Otherwise it's not symmetric.)
1202 31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	446 for i = dComps
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	447 Z = sparse(m_tot, m_tot);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	448 for l = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	449 for k = 1:dim
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	450 Z = Z + nu(l)C{l,i,k,j}nu(k);
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	451 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	452 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	453 Z = -tuningdim/h11Z;
1202 31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	454 X = Z;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	455 closure = closure + E{i}RHOiHiX'eH_gamma(e'*E{j}' );
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	456 penalty = penalty - E{i}RHOiHiX'e*H_gamma;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	457 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	458
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	459 % Free boundary condition
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	460 case {'F','f','Free','free','traction','Traction','t','T'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	461 closure = closure - E{j}RHOiHieH_gamma*tau';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	462 penalty = penalty + E{j}RHOiHieH_gamma;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	463
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	464 % Unknown boundary condition
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	465 otherwise
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	466 error('No such boundary condition: type = %s',type);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	467 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	468 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	469
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	470 % type Struct that specifies the interface coupling.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	471 % Fields:
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	472 % -- tuning: penalty strength, defaults to 1.0
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	473 % -- interpolation: type of interpolation, default 'none'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	474 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	475
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	476 defaultType.tuning = 1.0;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	477 defaultType.interpolation = 'none';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	478 default_struct('type', defaultType);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	479
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	480 switch type.interpolation
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	481 case {'none', ''}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	482 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	483 case {'op','OP'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	484 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	485 otherwise
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	486 error('Unknown type of interpolation: %s ', type.interpolation);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	487 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	488 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	489
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	490 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	491 tuning = type.tuning;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	492
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	493 % u denotes the solution in the own domain
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	494 % v denotes the solution in the neighbour domain
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	495
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	496 u = obj;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	497 v = neighbour_scheme;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	498
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	499 % Operators, u side
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	500 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	501 tau_u = u.getBoundaryOperator('tau', boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	502 h11_u = u.getBorrowing(boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	503 nu_u = u.getNormal(boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	504
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	505 E_u = u.E;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	506 C_u = u.C;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	507 m_tot_u = u.grid.N();
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	508
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	509 % Operators, v side
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	510 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	511 tau_v = v.getBoundaryOperator('tau', neighbour_boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	512 h11_v = v.getBorrowing(neighbour_boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	513 nu_v = v.getNormal(neighbour_boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	514
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	515 E_v = v.E;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	516 C_v = v.C;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	517 m_tot_v = v.grid.N();
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	518
1252 8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	519 % Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	520 flipFlag = false;
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	521 e_v_flip = e_v;
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	522 if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) \|\| ...
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	523 (strcmp(boundary,'e') && strcmp(neighbour_boundary,'s')) \|\| ...
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	524 (strcmp(boundary,'w') && strcmp(neighbour_boundary,'n')) \|\| ...
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	525 (strcmp(boundary,'n') && strcmp(neighbour_boundary,'w')) \|\| ...
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	526 (strcmp(boundary,'s') && strcmp(neighbour_boundary,'s')) \|\| ...
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	527 (strcmp(boundary,'n') && strcmp(neighbour_boundary,'n')) \|\| ...
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	528 (strcmp(boundary,'w') && strcmp(neighbour_boundary,'w')) \|\| ...
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	529 (strcmp(boundary,'e') && strcmp(neighbour_boundary,'e'))
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	530
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	531 flipFlag = true;
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	532 e_v_flip = fliplr(e_v);
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	533
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	534 t1 = tau_v(:,1:2:end-1);
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	535 t2 = tau_v(:,2:2:end);
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	536
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	537 t1 = fliplr(t1);
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	538 t2 = fliplr(t2);
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	539
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	540 tau_v(:,1:2:end-1) = t1;
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	541 tau_v(:,2:2:end) = t2;
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	542 end
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	543
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	544 % Operators that are only required for own domain
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	545 Hi = u.Hi_kron;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	546 RHOi = u.RHOi_kron;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	547 e_kron = u.getBoundaryOperator('e', boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	548 T_u = u.getBoundaryTractionOperator(boundary);
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	549
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	550 % Shared operators
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	551 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	552 H_gamma_kron = u.getBoundaryQuadrature(boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	553 dim = u.dim;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	554
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	555 % Preallocate
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	556 [~, m_int] = size(H_gamma);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	557 closure = sparse(dimm_tot_u, dimm_tot_u);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	558 penalty = sparse(dimm_tot_u, dimm_tot_v);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	559
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	560 % ---- Continuity of displacement ------
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	561
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	562 % Y: symmetrizing part of penalty
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	563 % Z: symmetric part of penalty
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	564 % X = Y + Z.
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	565
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	566 % Loop over components to couple across interface
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	567 for j = 1:dim
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	568
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	569 % Loop over components that penalties end up on
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	570 for i = 1:dim
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	571 Y = 1/2*T_u{j,i}';
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	572 Z_u = sparse(m_int, m_int);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	573 Z_v = sparse(m_int, m_int);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	574 for l = 1:dim
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	575 for k = 1:dim
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	576 Z_u = Z_u + e_u'nu_u(l)C_u{l,i,k,j}nu_u(k)e_u;
05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	577 Z_v = Z_v + e_v'nu_v(l)C_v{l,i,k,j}nu_v(k)e_v;
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	578 end
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	579 end
1252 8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	580
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	581 if flipFlag
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	582 Z_v = rot90(Z_v,2);
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	583 end
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	584
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	585 Z = -tuningdim( 1/(4h11_u)Z_u + 1/(4h11_v)Z_v );
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	586 X = Y + Z*e_u';
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	587 closure = closure + E_u{i}X'H_gammae_u'E_u{j}';
1252 8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	588 penalty = penalty - E_u{i}X'H_gammae_v_flip'E_v{j}';
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	589
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	590 end
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	591 end
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	592
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	593 % ---- Continuity of traction ------
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	594 closure = closure - 1/2e_kronH_gamma_kron*tau_u';
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	595 penalty = penalty - 1/2e_kronH_gamma_kron*tau_v';
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	596
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	597 % ---- Multiply by inverse of density x quadraure ----
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	598 closure = RHOiHiclosure;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	599 penalty = RHOiHipenalty;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	600
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	601 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	602
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	603 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	604 error('Non-conforming interfaces not implemented yet.');
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	605 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	606
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	607 % Returns the component number that is the tangential/normal component
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	608 % at the specified boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	609 function comp = getComponent(obj, comp_str, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	610 assertIsMember(comp_str, {'normal', 'tangential'});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	611 assertIsMember(boundary, {'w', 'e', 's', 'n'});
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	612
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	613 switch boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	614 case {'w', 'e'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	615 switch comp_str
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	616 case 'normal'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	617 comp = 1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	618 case 'tangential'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	619 comp = 2;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	620 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	621 case {'s', 'n'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	622 switch comp_str
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	623 case 'normal'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	624 comp = 2;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	625 case 'tangential'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	626 comp = 1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	627 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	628 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	629 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	630
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	631 % Returns h11 for the boundary specified by the string boundary.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	632 % op -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	633 function h11 = getBorrowing(obj, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	634 assertIsMember(boundary, {'w', 'e', 's', 'n'})
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	635
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	636 switch boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	637 case {'w','e'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	638 h11 = obj.h11{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	639 case {'s', 'n'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	640 h11 = obj.h11{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	641 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	642 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	643
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	644 % Returns the outward unit normal vector for the boundary specified by the string boundary.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	645 function nu = getNormal(obj, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	646 assertIsMember(boundary, {'w', 'e', 's', 'n'})
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	647
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	648 switch boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	649 case 'w'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	650 nu = [-1,0];
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	651 case 'e'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	652 nu = [1,0];
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	653 case 's'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	654 nu = [0,-1];
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	655 case 'n'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	656 nu = [0,1];
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	657 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	658 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	659
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	660 % Returns the boundary operator op for the boundary specified by the string boundary.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	661 % op -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	662 function o = getBoundaryOperator(obj, op, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	663 assertIsMember(boundary, {'w', 'e', 's', 'n'})
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	664 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'})
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	665
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	666 switch op
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	667 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	668 o = obj.([op, '_', boundary]);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	669 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	670
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	671 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	672
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	673 % Returns the boundary operator op for the boundary specified by the string boundary.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	674 % op -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	675 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	676 assertIsMember(boundary, {'w', 'e', 's', 'n'})
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	677 assertIsMember(op, {'e'})
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	678
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	679 switch op
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	680
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	681 case 'e'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	682 o = obj.(['e_scalar', '_', boundary]);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	683 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	684
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	685 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	686
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	687 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	688 % Formula: tau_i = T_ij u_j
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	689 % op -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	690 function T = getBoundaryTractionOperator(obj, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	691 assertIsMember(boundary, {'w', 'e', 's', 'n'})
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	692
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	693 T = obj.(['T', '_', boundary]);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	694 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	695
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	696 % Returns square boundary quadrature matrix, of dimension
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	697 % corresponding to the number of boundary unknowns
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	698 %
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	699 % boundary -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	700 function H = getBoundaryQuadrature(obj, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	701 assertIsMember(boundary, {'w', 'e', 's', 'n'})
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	702
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	703 H = obj.getBoundaryQuadratureForScalarField(boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	704 I_dim = speye(obj.dim, obj.dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	705 H = kron(H, I_dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	706 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	707
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	708 % Returns square boundary quadrature matrix, of dimension
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	709 % corresponding to the number of boundary grid points
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	710 %
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	711 % boundary -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	712 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	713 assertIsMember(boundary, {'w', 'e', 's', 'n'})
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	714
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	715 H_b = obj.(['H_', boundary]);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	716 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	717
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	718 function N = size(obj)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	719 N = obj.dim*prod(obj.m);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	720 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	721 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	722 end

Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dVariableAnisotropic.m @ 1253:89dad61cad22 feature/poroelastic