Mercurial > repos > public > sbplib

comparison +scheme/ViscoElastic2d.m @ 1308:5016f3f3badb feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add viscoElastic traction bc for normal-tangential
author Martin Almquist <malmquist@stanford.edu>
date Sun, 19 Jul 2020 21:24:48 -0700
parents fcca6ad8b102
children f59b849df30f
comparison
equal deleted inserted replaced
1307:fcca6ad8b102 1308:5016f3f3badb
38 38
39 % Bundary restriction ops 39 % Bundary restriction ops
40 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n 40 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n
41 41
42 n_w, n_e, n_s, n_n % Physical normals 42 n_w, n_e, n_s, n_n % Physical normals
43 tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents
43 44
44 elasticObj 45 elasticObj
45 46
46 end 47 end
47 48
205 206
206 obj.n_w = elasticObj.n_w; 207 obj.n_w = elasticObj.n_w;
207 obj.n_e = elasticObj.n_e; 208 obj.n_e = elasticObj.n_e;
208 obj.n_s = elasticObj.n_s; 209 obj.n_s = elasticObj.n_s;
209 obj.n_n = elasticObj.n_n; 210 obj.n_n = elasticObj.n_n;
211
212 obj.tangent_w = elasticObj.tangent_w;
213 obj.tangent_e = elasticObj.tangent_e;
214 obj.tangent_s = elasticObj.tangent_s;
215 obj.tangent_n = elasticObj.tangent_n;
210 216
211 obj.H_w = elasticObj.H_w; 217 obj.H_w = elasticObj.H_w;
212 obj.H_e = elasticObj.H_e; 218 obj.H_e = elasticObj.H_e;
213 obj.H_s = elasticObj.H_s; 219 obj.H_s = elasticObj.H_s;
214 obj.H_n = elasticObj.H_n; 220 obj.H_n = elasticObj.H_n;
265 eU = obj.eU; 271 eU = obj.eU;
266 eGamma = obj.eGamma; 272 eGamma = obj.eGamma;
267 273
268 switch type 274 switch type
269 case {'F','f','Free','free','traction','Traction','t','T'} 275 case {'F','f','Free','free','traction','Traction','t','T'}
276
277 % Get elastic closure and penalty
270 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning); 278 [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning);
271 closure = Eu*closure*Eu'; 279 closure = Eu*closure*Eu';
272 penalty = Eu*penalty; 280 penalty = Eu*penalty;
273 281
274 j = component; 282 switch component
275 for i = 1:dim 283 case 't'
276 for k = 1:dim 284 dir = obj.getTangent(boundary);
277 for l = 1:dim 285 case 'n'
278 closure = closure + eU{j}*( (RHO*H)\(C{i,j,k,l}*e*H_gamma*n{i}*e'*eGamma{k,l}') ); 286 dir = obj.getNormal(boundary);
279 end 287 case 1
280 end 288 dir = {1, 0};
281 end 289 case 2
290 dir = {0, 1};
291 end
292
293 % Add viscous part of closure
294 for j = 1:dim
295 for i = 1:dim
296 for k = 1:dim
297 for l = 1:dim
298 closure = closure + eU{j}*( (RHO*H)\(C{i,j,k,l}*e*dir{j}^2*H_gamma*n{i}*e'*eGamma{k,l}') );
299 end
300 end
301 end
302 end
303
282 end 304 end
283 305
284 end 306 end
285 307
286 function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning) 308 function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning)