sbplib: +scheme/ViscoElastic2d.m annotate

annotate +scheme/ViscoElastic2d.m @ 1314:58df4a35fe43 feature/poroelastic

Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.

author	Martin Almquist <malmquist@stanford.edu>
date	Sat, 25 Jul 2020 15:56:31 -0700
parents	a4c8bf2371f3
children	6650b111742d

rev	line source
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	1 classdef ViscoElastic2d < scheme.Scheme
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	2
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	3 % Discretizes the visco-elastic wave equation in curvilinear coordinates.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	4 % Assumes fully compatible operators.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	5
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	6 properties
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	7 m % Number of points in each direction, possibly a vector
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	8 h % Grid spacing
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	9
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	10 grid
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	11 dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	12
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	13 order % Order of accuracy for the approximation
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	14
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	15 % Diagonal matrices for variable coefficients
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	16 % J, Ji
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	17 RHO % Density
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	18 C % Elastic stiffness tensor
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	19 ETA % Effective viscosity, used in strain rate eq
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	20
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	21 D % Total operator
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	22 Delastic % Elastic operator (momentum balance)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	23 Dviscous % Acts on viscous strains in momentum balance
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	24 DstrainRate % Acts on u and gamma, returns strain rate gamma_t
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	25
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	26 D1, D1Tilde % Physical derivatives
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	27 sigma % Cell matrix of physical stress operators
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	28
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	29 % Inner products
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	30 H
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	31
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	32 % Boundary inner products (for scalar field)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	33 H_w, H_e, H_s, H_n
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	34
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	35 % Restriction operators
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	36 Eu, Egamma % Pick out all components of u/gamma
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	37 eU, eGamma % Pick out one specific component
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	38
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	39 % Bundary restriction ops
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	40 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	41
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	42 n_w, n_e, n_s, n_n % Physical normals
1308 5016f3f3badb Add viscoElastic traction bc for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1307 diff changeset	43 tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	44
1314 58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	45 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	46 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	47 tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	48 tau_t_w, tau_t_e, tau_t_s, tau_t_n % Return scalar field
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	49
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	50 elasticObj
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	51
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	52 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	53
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	54 methods
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	55
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	56 % The coefficients can either be function handles or grid functions
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	57 function obj = ViscoElastic2d(g, order, rho, C, eta)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	58 default_arg('rho', @(x,y) 0*x+1);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	59 default_arg('eta', @(x,y) 0*x+1);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	60 dim = 2;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	61
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	62 C_default = cell(dim,dim,dim,dim);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	63 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	64 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	65 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	66 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	67 C_default{i,j,k,l} = @(x,y) 0*x ;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	68 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	69 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	70 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	71 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	72 default_arg('C', C_default);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	73
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	74 assert(isa(g, 'grid.Curvilinear'));
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	75
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	76 if isa(rho, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	77 rho = grid.evalOn(g, rho);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	78 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	79
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	80 if isa(eta, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	81 eta = grid.evalOn(g, eta);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	82 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	83
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	84 C_mat = cell(dim,dim,dim,dim);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	85 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	86 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	87 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	88 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	89 if isa(C{i,j,k,l}, 'function_handle')
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	90 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	91 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	92 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	93 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	94 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	95 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	96 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	97 obj.C = C_mat;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	98
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	99 elasticObj = scheme.Elastic2dCurvilinearAnisotropic(g, order, rho, C);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	100
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	101 % Construct a pair of first derivatives
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	102 K = elasticObj.K;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	103 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	104 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	105 K{i,j} = spdiag(K{i,j});
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	106 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	107 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	108 J = elasticObj.J;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	109 Ji = elasticObj.Ji;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	110 D_ref = elasticObj.refObj.D1;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	111 D1 = cell(dim, 1);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	112 D1Tilde = cell(dim, 1);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	113 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	114 D1{i} = 0*D_ref{i};
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	115 D1Tilde{i} = 0*D_ref{i};
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	116 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	117 D1{i} = D1{i} + K{i,j}*D_ref{j};
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	118 D1Tilde{i} = D1Tilde{i} + JiD_ref{j}J*K{i,j};
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	119 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	120 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	121 obj.D1 = D1;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	122 obj.D1Tilde = D1Tilde;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	123
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	124 eU = elasticObj.E;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	125
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	126 % Storage order for gamma: 11-12-21-22.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	127 I = speye(g.N(), g.N());
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	128 eGamma = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	129 e = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	130 e{1,1} = [1;0;0;0];
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	131 e{1,2} = [0;1;0;0];
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	132 e{2,1} = [0;0;1;0];
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	133 e{2,2} = [0;0;0;1];
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	134 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	135 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	136 eGamma{i,j} = kron(I, e{i,j});
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	137 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	138 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	139
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	140 % Store u first, then gamma
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	141 mU = dim*g.N();
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	142 mGamma = dim^2*g.N();
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	143 Iu = speye(mU, mU);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	144 Igamma = speye(mGamma, mGamma);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	145
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	146 Eu = cell2mat({Iu, sparse(mU, mGamma)})';
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	147 Egamma = cell2mat({sparse(mGamma, mU), Igamma})';
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	148
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	149 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	150 eU{i} = Eu*eU{i};
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	151 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	152 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	153 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	154 eGamma{i,j} = Egamma*eGamma{i,j};
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	155 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	156 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	157
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	158 obj.eGamma = eGamma;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	159 obj.eU = eU;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	160 obj.Egamma = Egamma;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	161 obj.Eu = Eu;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	162
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	163 % Build stress operator
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	164 sigma = cell(dim, dim);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	165 C = obj.C;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	166 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	167 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	168 sigma{i,j} = spalloc(g.N(), (dim^2 + dim)g.N(), order^2g.N());
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	169 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	170 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	171 sigma{i,j} = sigma{i,j} + C{i,j,k,l}(D1{k}eU{l}' - eGamma{k,l}');
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	172 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	173 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	174 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	175 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	176
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	177 % Elastic operator
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	178 Delastic = EuelasticObj.DEu';
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	179
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	180 % Add viscous strains to momentum balance
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	181 RHOi = spdiag(1./rho);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	182 Dviscous = spalloc((dim^2 + dim)g.N(), (dim^2 + dim)g.N(), order^2(dim^2 + dim)g.N());
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	183 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	184 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	185 for k = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	186 for l = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	187 Dviscous = Dviscous - eU{j}RHOiD1Tilde{i}C{i,j,k,l}eGamma{k,l}';
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	188 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	189 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	190 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	191 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	192
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	193 ETA = spdiag(eta);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	194 DstrainRate = 0*Delastic;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	195 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	196 for j = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	197 DstrainRate = DstrainRate + eGamma{i,j}*(ETA\sigma{i,j});
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	198 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	199 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	200
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	201 obj.D = Delastic + Dviscous + DstrainRate;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	202 obj.Delastic = Delastic;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	203 obj.Dviscous = Dviscous;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	204 obj.DstrainRate = DstrainRate;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	205 obj.sigma = sigma;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	206
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	207 %---- Set remaining object properties ------
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	208 obj.RHO = elasticObj.RHO;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	209 obj.ETA = ETA;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	210 obj.H = elasticObj.H;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	211
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	212 obj.n_w = elasticObj.n_w;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	213 obj.n_e = elasticObj.n_e;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	214 obj.n_s = elasticObj.n_s;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	215 obj.n_n = elasticObj.n_n;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	216
1308 5016f3f3badb Add viscoElastic traction bc for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1307 diff changeset	217 obj.tangent_w = elasticObj.tangent_w;
5016f3f3badb Add viscoElastic traction bc for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1307 diff changeset	218 obj.tangent_e = elasticObj.tangent_e;
5016f3f3badb Add viscoElastic traction bc for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1307 diff changeset	219 obj.tangent_s = elasticObj.tangent_s;
5016f3f3badb Add viscoElastic traction bc for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1307 diff changeset	220 obj.tangent_n = elasticObj.tangent_n;
5016f3f3badb Add viscoElastic traction bc for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1307 diff changeset	221
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	222 obj.H_w = elasticObj.H_w;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	223 obj.H_e = elasticObj.H_e;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	224 obj.H_s = elasticObj.H_s;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	225 obj.H_n = elasticObj.H_n;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	226
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	227 obj.e_scalar_w = elasticObj.e_scalar_w;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	228 obj.e_scalar_e = elasticObj.e_scalar_e;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	229 obj.e_scalar_s = elasticObj.e_scalar_s;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	230 obj.e_scalar_n = elasticObj.e_scalar_n;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	231
1314 58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	232 % -- Create new traction operators including viscous strain contribution --
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	233 tau1 = struct;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	234 tau2 = struct;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	235 tau_n = struct;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	236 tau_t = struct;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	237 boundaries = {'w', 'e', 's', 'n'};
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	238 for bNumber = 1:numel(boundaries)
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	239 b = boundaries{bNumber};
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	240
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	241 n = elasticObj.getNormal(b);
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	242 t = elasticObj.getTangent(b);
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	243 e = elasticObj.getBoundaryOperatorForScalarField('e', b);
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	244 tau1.(b) = (elasticObj.getBoundaryOperator('tau1', b)'*Eu')';
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	245 tau2.(b) = (elasticObj.getBoundaryOperator('tau2', b)'*Eu')';
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	246
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	247 % Add viscous contributions
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	248 for i = 1:dim
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	249 for k = 1:dim
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	250 for l = 1:dim
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	251 tau1.(b) = tau1.(b) - (n{i}e'C{i,1,k,l}*eGamma{k,l}')';
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	252 tau2.(b) = tau2.(b) - (n{i}e'C{i,2,k,l}*eGamma{k,l}')';
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	253 end
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	254 end
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	255 end
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	256
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	257 % Compute normal and tangential components
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	258 tau_n.(b) = tau1.(b)n{1} + tau2.(b)n{2};
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	259 tau_t.(b) = tau1.(b)t{1} + tau2.(b)t{2};
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	260 end
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	261
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	262 obj.tau1_w = tau1.w;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	263 obj.tau1_e = tau1.e;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	264 obj.tau1_s = tau1.s;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	265 obj.tau1_n = tau1.n;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	266
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	267 obj.tau2_w = tau2.w;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	268 obj.tau2_e = tau2.e;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	269 obj.tau2_s = tau2.s;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	270 obj.tau2_n = tau2.n;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	271
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	272 obj.tau_n_w = tau_n.w;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	273 obj.tau_n_e = tau_n.e;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	274 obj.tau_n_s = tau_n.s;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	275 obj.tau_n_n = tau_n.n;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	276
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	277 obj.tau_t_w = tau_t.w;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	278 obj.tau_t_e = tau_t.e;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	279 obj.tau_t_s = tau_t.s;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	280 obj.tau_t_n = tau_t.n;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	281 %----------------------------------------
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	282
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	283 % Misc.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	284 obj.elasticObj = elasticObj;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	285 obj.m = elasticObj.m;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	286 obj.h = elasticObj.h;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	287
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	288 obj.order = order;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	289 obj.grid = g;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	290 obj.dim = dim;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	291
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	292 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	293
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	294
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	295 % Closure functions return the operators applied to the own domain to close the boundary
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	296 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	297 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	298 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	299 % on the first component. Can also be e.g.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	300 % {'normal', 'd'} or {'tangential', 't'} for conditions on
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	301 % tangential/normal component.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	302 % data is a function returning the data that should be applied at the boundary.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	303 % neighbour_scheme is an instance of Scheme that should be interfaced to.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	304 % neighbour_boundary is a string specifying which boundary to interface to.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	305
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	306 % For displacement bc:
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	307 % bc = {comp, 'd', dComps},
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	308 % where
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	309 % dComps = vector of components with displacement BC. Default: 1:dim.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	310 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	311 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	312 default_arg('tuning', 1.0);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	313 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	314
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	315 component = bc{1};
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	316 type = bc{2};
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	317 dim = obj.dim;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	318
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	319 n = obj.getNormal(boundary);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	320 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	321 e = obj.getBoundaryOperatorForScalarField('e', boundary);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	322
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	323 H = obj.H;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	324 RHO = obj.RHO;
1309 f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	325 ETA = obj.ETA;
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	326 C = obj.C;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	327 Eu = obj.Eu;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	328 eU = obj.eU;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	329 eGamma = obj.eGamma;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	330
1311 b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	331 % Get elastic closure and penalty
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	332 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning);
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	333 closure = EuclosureEu';
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	334 penalty = Eu*penalty;
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	335
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	336 switch component
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	337 case 't'
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	338 dir = obj.getTangent(boundary);
1314 58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	339 tau = obj.getBoundaryOperator('tau_t', boundary);
1311 b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	340 case 'n'
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	341 dir = obj.getNormal(boundary);
1314 58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	342 tau = obj.getBoundaryOperator('tau_n', boundary);
1311 b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	343 case 1
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	344 dir = {1, 0};
1314 58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	345 tau = obj.getBoundaryOperator('tau1', boundary);
1311 b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	346 case 2
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	347 dir = {0, 1};
1314 58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	348 tau = obj.getBoundaryOperator('tau2', boundary);
1311 b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	349 end
b2e84fe336d8 Move common code for displacement and traction BC outside of switch statement Martin Almquist <malmquist@stanford.edu> parents: 1310 diff changeset	350
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	351 switch type
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	352 case {'F','f','Free','free','traction','Traction','t','T'}
1308 5016f3f3badb Add viscoElastic traction bc for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1307 diff changeset	353
1314 58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	354 % Set elastic closure to zero
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	355 closure = 0*closure;
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	356
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	357 for m = 1:dim
58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	358 closure = closure - eU{m}( (RHOH)\(edir{m}H_gamma*tau') );
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	359 end
1308 5016f3f3badb Add viscoElastic traction bc for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1307 diff changeset	360
1309 f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	361 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	362
f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	363 % Add penalty to strain rate eq
f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	364 for i = 1:dim
f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	365 for j = 1:dim
f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	366 for k = 1:dim
1310 eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1309 diff changeset	367 for l = 1:dim
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1309 diff changeset	368 for m = 1:dim
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1309 diff changeset	369 closure = closure - eGamma{i,j}( (HETA)\(C{i,j,k,l}eH_gammadir{l}dir{m}n{k}e'*eU{m}') );
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1309 diff changeset	370 end
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1309 diff changeset	371 penalty = penalty + eGamma{i,j}( (HETA)\(C{i,j,k,l}eH_gammadir{l}n{k}) );
eb015fe9605b Add viscoelastic displacement BC for normal-tangential. Bugfix traction BC for normal-tangential Martin Almquist <malmquist@stanford.edu> parents: 1309 diff changeset	372 end
1309 f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	373 end
f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	374 end
f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	375 end
f59b849df30f Add displacement BC to viscoelastic Martin Almquist <malmquist@stanford.edu> parents: 1308 diff changeset	376
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	377 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	378
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	379 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	380
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	381 % type Struct that specifies the interface coupling.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	382 % Fields:
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	383 % -- tuning: penalty strength, defaults to 1.0
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	384 % -- interpolation: type of interpolation, default 'none'
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	385 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	386
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	387 defaultType.tuning = 1.0;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	388 defaultType.interpolation = 'none';
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	389 defaultType.type = 'standard';
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	390 default_struct('type', defaultType);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	391
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	392 switch type.type
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	393 case 'standard'
1312 6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	394 [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type);
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	395 case 'frictionalFault'
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	396 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	397 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	398
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	399 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	400
1312 6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	401 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	402
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	403 % u denotes the solution in the own domain
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	404 % v denotes the solution in the neighbour domain
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	405 u = obj;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	406 v = neighbour_scheme;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	407
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	408 dim = obj.dim;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	409
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	410 n = u.getNormal(boundary);
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	411 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	412 e = u.getBoundaryOperatorForScalarField('e', boundary);
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	413
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	414 ev = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	415
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	416 H = u.H;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	417 RHO = u.RHO;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	418 ETA = u.ETA;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	419 C = u.C;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	420 Eu = u.Eu;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	421 eU = u.eU;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	422 eGamma = u.eGamma;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	423
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	424 CV = v.C;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	425 Ev = v.Eu;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	426 eV = v.eU;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	427 eGammaV = v.eGamma;
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	428 nV = v.getNormal(neighbour_boundary);
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	429
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	430
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	431 % Get elastic closure and penalty
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	432 [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	433 closure = EuclosureEu';
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	434 penalty = EupenaltyEv';
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	435
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	436 % Add viscous part of traction coupling
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	437 for i = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	438 for j = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	439 for k = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	440 for l = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	441 closure = closure + 1/2eU{j}( (RHOH)\(C{i,j,k,l}eH_gamman{i}e'eGamma{k,l}') );
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	442 penalty = penalty + 1/2eU{j}( (RHOH)\(eH_gammanV{i}(ev'CV{i,j,k,l}ev)ev'eGammaV{k,l}') );
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	443 end
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	444 end
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	445 end
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	446 end
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	447
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	448 % Add penalty to strain rate eq
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	449 for i = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	450 for j = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	451 for k = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	452 for l = 1:dim
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	453 closure = closure - 1/2eGamma{i,j}( (HETA)\(C{i,j,k,l}eH_gamman{k}e'eU{l}') );
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	454 penalty = penalty + 1/2eGamma{i,j}( (HETA)\(C{i,j,k,l}eH_gamman{k}ev'eV{l}') );
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	455 end
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	456 end
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	457 end
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	458 end
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	459
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	460
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	461 end
6d64b57caf46 Add viscoelastic regular interfaces Martin Almquist <malmquist@stanford.edu> parents: 1311 diff changeset	462
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	463 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	464 tuning = type.tuning;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	465
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	466 % u denotes the solution in the own domain
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	467 % v denotes the solution in the neighbour domain
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	468 u = obj;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	469 v = neighbour_scheme;
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	470
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	471 dim = obj.dim;
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	472
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	473 n = u.getNormal(boundary);
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	474 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	475 e = u.getBoundaryOperatorForScalarField('e', boundary);
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	476
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	477 ev = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	478
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	479 H = u.H;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	480 RHO = u.RHO;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	481 ETA = u.ETA;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	482 C = u.C;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	483 Eu = u.Eu;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	484 eU = u.eU;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	485 eGamma = u.eGamma;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	486 Egamma = u.Egamma;
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	487
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	488 CV = v.C;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	489 Ev = v.Eu;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	490 eV = v.eU;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	491 eGammaV = v.eGamma;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	492 nV = v.getNormal(neighbour_boundary);
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	493
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	494 % Reduce stiffness tensors to boundary size
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	495 for i = 1:dim
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	496 for j = 1:dim
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	497 for k = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	498 for l = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	499 C{i,j,k,l} = e'C{i,j,k,l}e;
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	500 CV{i,j,k,l} = ev'CV{i,j,k,l}ev;
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	501 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	502 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	503 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	504 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	505
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	506 % Get elastic closure and penalty
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	507 [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	508 closure = EuclosureEu';
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	509 penalty = EupenaltyEv';
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	510
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	511 % ---- Tangential tractions are imposed just like traction BC ------
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	512 % We only need the viscous part
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	513 closure_tangential = obj.boundary_condition(boundary, {'t', 't'});
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	514 closure = closure + closure_tangentialEgammaEgamma';
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	515
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	516
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	517 % ------ Coupling of normal component -----------
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	518 % Add viscous part of traction coupling
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	519 for i = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	520 for j = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	521 for k = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	522 for l = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	523 for m = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	524 closure = closure + 1/2eU{m}( (RHOH)\(en{m}H_gamman{j}n{i}C{i,j,k,l}e'eGamma{k,l}') );
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	525 penalty = penalty - 1/2eU{m}( (RHOH)\(en{m}H_gammanV{j}nV{i}CV{i,j,k,l}ev'eGammaV{k,l}') );
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	526 end
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	527 end
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	528 end
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	529 end
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	530 end
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	531
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	532 % Add penalty to strain rate eq
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	533 for i = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	534 for j = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	535 for k = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	536 for l = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	537 for m = 1:dim
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	538 closure = closure - 1/2eGamma{i,j}( (HETA)\(en{l}n{k}C{i,j,k,l}H_gamman{m}e'eU{m}') );
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	539 penalty = penalty - 1/2eGamma{i,j}( (HETA)\(en{l}n{k}C{i,j,k,l}H_gammanV{m}ev'eV{m}') );
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	540 end
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	541 end
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	542 end
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	543 end
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	544 end
1313 a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	545 %-------------------------------------------------
a4c8bf2371f3 Add viscoelastic fault interface Martin Almquist <malmquist@stanford.edu> parents: 1312 diff changeset	546
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	547 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	548
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	549 % Returns the outward unit normal vector for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	550 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	551 function n = getNormal(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	552 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	553
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	554 n = obj.(['n_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	555 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	556
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	557 % Returns the unit tangent vector for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	558 % t is a cell of diagonal matrices for each normal component, t{1} = t_1, t{2} = t_2.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	559 function t = getTangent(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	560 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	561
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	562 t = obj.(['tangent_' boundary]);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	563 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	564
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	565 % Returns the boundary operator op for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	566 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	567 function o = getBoundaryOperator(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	568 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	569 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	570
1314 58df4a35fe43 Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC. Martin Almquist <malmquist@stanford.edu> parents: 1313 diff changeset	571 o = obj.([op, '_', boundary]);
1307 fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	572
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	573 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	574
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	575 % Returns the boundary operator op for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	576 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	577 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	578 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	579 assertIsMember(op, {'e'})
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	580
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	581 switch op
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	582
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	583 case 'e'
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	584 o = obj.(['e_scalar', '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	585 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	586
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	587 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	588
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	589 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	590 % Formula: tau_i = T_ij u_j
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	591 % op -- string
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	592 function T = getBoundaryTractionOperator(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	593 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	594
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	595 T = obj.(['T', '_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	596 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	597
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	598 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	599 % corresponding to the number of boundary unknowns
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	600 %
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	601 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	602 function H = getBoundaryQuadrature(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	603 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	604
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	605 H = obj.getBoundaryQuadratureForScalarField(boundary);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	606 I_dim = speye(obj.dim, obj.dim);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	607 H = kron(H, I_dim);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	608 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	609
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	610 % Returns square boundary quadrature matrix, of dimension
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	611 % corresponding to the number of boundary grid points
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	612 %
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	613 % boundary -- string
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	614 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	615 assertIsMember(boundary, {'w', 'e', 's', 'n'})
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	616
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	617 H_b = obj.(['H_', boundary]);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	618 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	619
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	620 function N = size(obj)
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	621 N = (obj.dim + obj.dim^2)*prod(obj.m);
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	622 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	623 end
fcca6ad8b102 Add diffOp for viscoElastic Martin Almquist <malmquist@stanford.edu> parents: diff changeset	624 end

Mercurial > repos > public > sbplib

annotate +scheme/ViscoElastic2d.m @ 1314:58df4a35fe43 feature/poroelastic