sbplib: +scheme/Elastic2dVariableAnisotropic.m annotate

annotate +scheme/Elastic2dVariableAnisotropic.m @ 1305:b5907140c069 feature/poroelastic

Clean up hollow D2 generation by using precomputed coefficients stored in mat files

author	Martin Almquist <malmquist@stanford.edu>
date	Mon, 13 Jul 2020 14:21:41 -0700
parents	a38e80fdbf60
children	8d9fc7981796

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.
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	58 function obj = Elastic2dVariableAnisotropic(g, order, rho, C, opSet, optFlag, hollow)
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	59 default_arg('hollow', false);
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	60 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	61 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	62 default_arg('optFlag', false);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	63 dim = 2;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	64
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	65 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	66 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	67 for j = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	68 for k = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	69 for l = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	70 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	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 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	75 default_arg('C', C_default);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	76 assert(isa(g, 'grid.Cartesian'))
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	77
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	78 if isa(rho, 'function_handle')
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	79 rho = grid.evalOn(g, rho);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	80 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	81
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	82 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	83 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	84 for j = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	85 for k = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	86 for l = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	87 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	88 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	89 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	90 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	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 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	95 obj.C = C_mat;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	96
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	97 m = g.size();
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	98 m_tot = g.N();
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	99 lim = g.lim;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	100 if isempty(lim)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	101 x = g.x;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	102 lim = cell(length(x),1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	103 for i = 1:length(x)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	104 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	105 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	106 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	107
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	108 % 1D operators
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	109 ops = cell(dim,1);
1303 49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	110 opsHollow = cell(dim,1);
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	111 h = zeros(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	112 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	113 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	114 h(i) = ops{i}.h;
1303 49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	115 if hollow
49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	116 opsHollow{i} = sbp.D2VariableCompatibleHollow(m(i), lim{i}, order);
49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	117 end
1201 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 % Borrowing constants
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	121 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	122 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	123 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	124
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	125 I = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	126 D1 = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	127 D2 = cell(dim,1);
1303 49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	128 D2Hollow = cell(dim,1);
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	129 H = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	130 Hi = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	131 e_0 = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	132 e_m = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	133 d1_0 = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	134 d1_m = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	135
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	136 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	137 I{i} = speye(m(i));
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	138 D1{i} = ops{i}.D1;
1303 49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	139 if hollow
49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	140 D2Hollow{i} = opsHollow{i}.D2;
49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	141 end
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	142 D2{i} = ops{i}.D2;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	143 H{i} = ops{i}.H;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	144 Hi{i} = ops{i}.HI;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	145 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	146 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	147 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	148 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	149 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	150
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	151 %====== Assemble full operators ========
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	152 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	153 RHO = spdiag(rho);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	154 obj.RHO = RHO;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	155 obj.RHOi = inv(RHO);
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	156 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	157
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	158 obj.D1 = cell(dim,1);
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	159 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	160
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	161 % D1
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	162 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	163 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	164
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	165 % Boundary restriction operators
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	166 e_l = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	167 e_r = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	168 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	169 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	170 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	171 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	172
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	173 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	174 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	175 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	176 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	177
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	178 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	179 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	180 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	181 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	182
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	183 % 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	184 E = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	185 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	186 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	187 e = sparse(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	188 e(i) = 1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	189 E{i} = kron(I,e);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	190 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	191 obj.E = E;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	192
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	193 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	194 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	195 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	196 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	197
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	198 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	199 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	200 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	201 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	202
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	203
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	204 % D2
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	205 switch order
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	206 case 2
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	207 width = 3;
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	208 nBP = 2;
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	209 case 4
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	210 width = 5;
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	211 nBP = 6;
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	212 case 6
89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	213 width = 7;
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	214 nBP = 9;
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	215 end
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	216 for j = 1:dim
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	217 for k = 1:dim
6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	218 for l = 1:dim
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	219 if hollow
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	220 D2_temp{j,k,l} = sparse(m_tot, m_tot);
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	221 else
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	222 D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot);
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	223 end
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	224 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	225 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	226 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	227 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	228
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	229 k = 1;
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	230 if hollow
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	231 mask = sparse(m(1), m(1));
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	232 mask(1:nBP, 1:nBP) = speye(nBP, nBP);
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	233 mask(end-nBP+1:end, end-nBP+1:end) = speye(nBP, nBP);
1304 a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	234 maskXSmall = kron(mask, speye(m(2), m(2)));
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	235 maskX = E{1}maskXSmallE{1}' + E{2}maskXSmallE{2}';
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	236 end
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	237 for r = 1:m(2)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	238 p = ind(:,r);
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	239 for j = 1:dim
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	240 for l = 1:dim
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	241 coeff = C{k,j,k,l};
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	242 if hollow && r > nBP && r < m(2) - nBP + 1
1303 49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	243 D_kk = D2Hollow{1}(coeff(p));
49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	244 else
49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	245 D_kk = D2{1}(coeff(p));
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	246 end
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	247 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	248 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	249 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	250 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	251
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	252 k = 2;
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	253 if hollow
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	254 mask = sparse(m(2), m(2));
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	255 mask(1:nBP, 1:nBP) = speye(nBP, nBP);
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	256 mask(end-nBP+1:end, end-nBP+1:end) = speye(nBP, nBP);
1304 a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	257 maskYSmall = kron(speye(m(1), m(1)), mask);
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	258
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	259 maskY = E{1}maskYSmallE{1}' + E{2}maskYSmallE{2}';
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	260 mask = maskX + maskY;
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	261 mask = mask>0;
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	262
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	263 maskSmall = maskXSmall + maskYSmall;
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	264 maskSmall = maskSmall>0;
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	265 end
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	266 for r = 1:m(1)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	267 p = ind(r,:);
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	268 for j = 1:dim
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	269 for l = 1:dim
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	270 coeff = C{k,j,k,l};
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	271 if hollow && r > nBP && r < m(1) - nBP + 1
1303 49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	272 D_kk = D2Hollow{2}(coeff(p));
49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	273 else
49e3870335ef Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet Martin Almquist <malmquist@stanford.edu> parents: 1302 diff changeset	274 D_kk = D2{2}(coeff(p));
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	275 end
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	276 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	277 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	278 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	279 end
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 % Quadratures
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	282 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	283 obj.Hi = inv(obj.H);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	284 obj.H_w = H{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	285 obj.H_e = H{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	286 obj.H_s = H{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	287 obj.H_n = H{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	288 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	289
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	290 % Differentiation matrix D (without SAT)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	291 D1 = obj.D1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	292 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	293 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	294 for j = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	295 for k = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	296 for l = 1:dim
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	297 if i == k
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	298 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	299 D2_temp{j,k,l} = [];
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	300 else
1304 a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	301 if hollow
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	302 D = D + E{j}(maskSmallD1{i})C_mat{i,j,k,l}D1{k}*E{l}';
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	303 else
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	304 D = D + E{j}(D1{i})C_mat{i,j,k,l}D1{k}E{l}';
a38e80fdbf60 Improve performance for the hollow case Martin Almquist <malmquist@stanford.edu> parents: 1303 diff changeset	305 end
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	306 end
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 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	309 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	310 end
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	311 clear D2_temp;
1302 a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	312 if hollow
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	313 mask = maskX + maskY;
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	314 mask = mask>0;
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	315 D = mask*D;
a0d615bde7f8 Add the hollow option to the anisotropic diffops Martin Almquist <malmquist@stanford.edu> parents: 1260 diff changeset	316 end
1206 6b203030fb37 Fix performance (full matrices) issues in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	317 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	318 obj.D = D;
1253 89dad61cad22 Make Elastic2dVariable faster and more memory efficient Martin Almquist <malmquist@stanford.edu> parents: 1252 diff changeset	319 clear D;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	320 %=========================================%'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	321
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	322 % Numerical traction operators for BC.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	323 %
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	324 % 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	325 %
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	326 T_l = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	327 T_r = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	328 tau_l = cell(dim,1);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	329 tau_r = cell(dim,1);
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 D1 = obj.D1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	332
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	333 % Boundary j
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	334 for j = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	335 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	336 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	337 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	338 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	339
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	340 [~, 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	341 [~, 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	342
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	343 % Traction component i
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	344 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	345 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	346 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	347
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	348 % Displacement component l
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	349 for l = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	350 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	351 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	352
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	353 % Derivative direction k
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	354 for k = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	355 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	356 - (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	357 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	358 + (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	359 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	360 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	361 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	362 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	363 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	364 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	365
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	366 % Traction tensors, T_ij
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	367 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	368 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	369 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	370 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	371
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	372 % Restriction operators
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	373 obj.e_w = e_w;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	374 obj.e_e = e_e;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	375 obj.e_s = e_s;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	376 obj.e_n = e_n;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	377
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	378 obj.e1_w = e1_w;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	379 obj.e1_e = e1_e;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	380 obj.e1_s = e1_s;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	381 obj.e1_n = e1_n;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	382
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	383 obj.e2_w = e2_w;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	384 obj.e2_e = e2_e;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	385 obj.e2_s = e2_s;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	386 obj.e2_n = e2_n;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	387
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	388 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	389 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	390 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	391 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	392
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	393 % First component of traction
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	394 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	395 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	396 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	397 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	398
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	399 % Second component of traction
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	400 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	401 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	402 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	403 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	404
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	405 % Traction vectors
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	406 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	407 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	408 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	409 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	410
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	411 % Kroneckered norms and coefficients
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	412 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	413
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	414 % Misc.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	415 obj.m = m;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	416 obj.h = h;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	417 obj.order = order;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	418 obj.grid = g;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	419 obj.dim = dim;
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 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	422
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	423
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	424 % 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	425 % 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	426 % 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	427 % 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	428 % 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	429 % {'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	430 % tangential/normal component.
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	431 % 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	432 % 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	433 % 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	434
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	435 % For displacement bc:
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	436 % bc = {comp, 'd', dComps},
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	437 % where
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	438 % dComps = vector of components with displacement BC. Default: 1:dim.
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	439 % 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	440 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	441 default_arg('tuning', 1.0);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	442
1202 31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	443 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	444 comp = bc{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	445 type = bc{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	446 if ischar(comp)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	447 comp = obj.getComponent(comp, boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	448 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	449
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	450 e = obj.getBoundaryOperatorForScalarField('e', boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	451 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	452 T = obj.getBoundaryTractionOperator(boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	453 h11 = obj.getBorrowing(boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	454 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	455 nu = obj.getNormal(boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	456
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	457 E = obj.E;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	458 Hi = obj.Hi;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	459 RHOi = obj.RHOi;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	460 C = obj.C;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	461
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	462 dim = obj.dim;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	463 m_tot = obj.grid.N();
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	464
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	465 % Preallocate
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	466 [~, col] = size(tau);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	467 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	468 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	469
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	470 j = comp;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	471 switch type
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	472
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	473 % Dirichlet boundary condition
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	474 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	475
1202 31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	476 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	477 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	478 else
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	479 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	480 end
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	481
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	482 % 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	483 % Y: symmetrizing part of penalty
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	484 % Z: symmetric part of penalty
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	485 % 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	486
31d7288d0653 Make Anisotropic work for mixed displacement/traction BC at the same boundary. Martin Almquist <malmquist@stanford.edu> parents: 1201 diff changeset	487 % Nonsymmetric part goes on all components to
1203 25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	488 % 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	489 for i = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	490 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	491 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	492 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	493 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	494 end
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	495
1203 25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	496 % Symmetric part only required on components with displacement BC.
25cadc69a589 Update comments in Anisotropic Martin Almquist <malmquist@stanford.edu> parents: 1202 diff changeset	497 % (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	498 for i = dComps
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	499 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	500 for l = 1:dim
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	501 for k = 1:dim
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	502 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	503 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	504 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	505 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	506 X = Z;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	507 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	508 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	509 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	510
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	511 % Free boundary condition
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	512 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	513 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	514 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	515
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	516 % Unknown boundary condition
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	517 otherwise
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	518 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	519 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	520 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	521
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	522 % 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	523 % Fields:
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	524 % -- 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	525 % -- 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	526 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	527
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	528 defaultType.tuning = 1.0;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	529 defaultType.interpolation = 'none';
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	530 default_struct('type', defaultType);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	531
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	532 switch type.interpolation
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	533 case {'none', ''}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	534 [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	535 case {'op','OP'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	536 [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	537 otherwise
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	538 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	539 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	540 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	541
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	542 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	543 tuning = type.tuning;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	544
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	545 % 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	546 % v denotes the solution in the neighbour domain
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	547
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	548 u = obj;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	549 v = neighbour_scheme;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	550
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	551 % Operators, u side
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	552 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	553 tau_u = u.getBoundaryOperator('tau', boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	554 h11_u = u.getBorrowing(boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	555 nu_u = u.getNormal(boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	556
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	557 E_u = u.E;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	558 C_u = u.C;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	559 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	560
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	561 % Operators, v side
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	562 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	563 tau_v = v.getBoundaryOperator('tau', neighbour_boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	564 h11_v = v.getBorrowing(neighbour_boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	565 nu_v = v.getNormal(neighbour_boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	566
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	567 E_v = v.E;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	568 C_v = v.C;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	569 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	570
1252 8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	571 % 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	572 flipFlag = false;
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	573 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	574 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	575 (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	576 (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	577 (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	578 (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	579 (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	580 (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	581 (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	582
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	583 flipFlag = true;
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	584 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	585
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	586 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	587 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	588
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	589 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	590 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	591
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	592 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	593 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	594 end
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	595
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	596 % Operators that are only required for own domain
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	597 Hi = u.Hi_kron;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	598 RHOi = u.RHOi_kron;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	599 e_kron = u.getBoundaryOperator('e', boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	600 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	601
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	602 % Shared operators
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	603 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	604 H_gamma_kron = u.getBoundaryQuadrature(boundary);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	605 dim = u.dim;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	606
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	607 % Preallocate
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	608 [~, m_int] = size(H_gamma);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	609 closure = sparse(dimm_tot_u, dimm_tot_u);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	610 penalty = sparse(dimm_tot_u, dimm_tot_v);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	611
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	612 % ---- Continuity of displacement ------
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	613
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	614 % Y: symmetrizing part of penalty
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	615 % Z: symmetric part of penalty
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	616 % X = Y + Z.
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	617
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	618 % Loop over components to couple across interface
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	619 for j = 1:dim
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	620
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	621 % Loop over components that penalties end up on
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	622 for i = 1:dim
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	623 Y = 1/2*T_u{j,i}';
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	624 Z_u = sparse(m_int, m_int);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	625 Z_v = sparse(m_int, m_int);
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	626 for l = 1:dim
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	627 for k = 1:dim
1207 05a01f77d0e3 Swap indices and fix bug in ElasticVariableAnisotropic Martin Almquist <malmquist@stanford.edu> parents: 1204 diff changeset	628 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	629 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	630 end
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	631 end
1252 8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	632
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	633 if flipFlag
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	634 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	635 end
8fc2f9a4c882 Add (temporary) fix for e-s , w-n, and x-x couplings Martin Almquist <malmquist@stanford.edu> parents: 1208 diff changeset	636
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	637 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	638 X = Y + Z*e_u';
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	639 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	640 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	641
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	642 end
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	643 end
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	644
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	645 % ---- Continuity of traction ------
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	646 closure = closure - 1/2e_kronH_gamma_kron*tau_u';
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	647 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	648
1204 687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	649 % ---- Multiply by inverse of density x quadraure ----
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	650 closure = RHOiHiclosure;
687515778437 Add anisotropic interface Martin Almquist <malmquist@stanford.edu> parents: 1203 diff changeset	651 penalty = RHOiHipenalty;
1201 8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	652
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	653 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	654
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	655 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	656 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	657 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	658
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	659 % 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	660 % at the specified boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	661 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	662 assertIsMember(comp_str, {'normal', 'tangential'});
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
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	665 switch boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	666 case {'w', 'e'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	667 switch comp_str
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	668 case 'normal'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	669 comp = 1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	670 case 'tangential'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	671 comp = 2;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	672 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	673 case {'s', 'n'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	674 switch comp_str
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	675 case 'normal'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	676 comp = 2;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	677 case 'tangential'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	678 comp = 1;
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	679 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	680 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	681 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	682
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	683 % 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	684 % op -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	685 function h11 = getBorrowing(obj, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	686 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	687
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	688 switch boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	689 case {'w','e'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	690 h11 = obj.h11{1};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	691 case {'s', 'n'}
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	692 h11 = obj.h11{2};
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	693 end
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 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	697 function nu = getNormal(obj, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	698 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	699
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	700 switch boundary
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	701 case 'w'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	702 nu = [-1,0];
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	703 case 'e'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	704 nu = [1,0];
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	705 case 's'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	706 nu = [0,-1];
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	707 case 'n'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	708 nu = [0,1];
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	709 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	710 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	711
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	712 % 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	713 % op -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	714 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	715 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	716 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	717
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	718 switch op
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	719 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	720 o = obj.([op, '_', boundary]);
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
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	723 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	724
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	725 % 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	726 % op -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	727 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	728 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	729 assertIsMember(op, {'e'})
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	730
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	731 switch op
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	732
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	733 case 'e'
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	734 o = obj.(['e_scalar', '_', boundary]);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	735 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	736
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	737 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	738
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	739 % 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	740 % 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	741 % op -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	742 function T = getBoundaryTractionOperator(obj, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	743 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	744
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	745 T = obj.(['T', '_', boundary]);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	746 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	747
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	748 % 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	749 % 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	750 %
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	751 % boundary -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	752 function H = getBoundaryQuadrature(obj, boundary)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	753 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	754
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	755 H = obj.getBoundaryQuadratureForScalarField(boundary);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	756 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	757 H = kron(H, I_dim);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	758 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	759
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	760 % 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	761 % 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	762 %
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	763 % boundary -- string
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	764 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	765 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	766
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	767 H_b = obj.(['H_', boundary]);
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	768 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	769
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	770 function N = size(obj)
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	771 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	772 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	773 end
8f4e79aa32ba Add fully compatible D2Variable opSet and first implementation of ElasticAnisotropic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	774 end

Mercurial > repos > public > sbplib

annotate +scheme/Elastic2dVariableAnisotropic.m @ 1305:b5907140c069 feature/poroelastic