Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropic.m @ 1291:a8e730db76e9 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add normal-tangential traction BC ElasticCurvilinearAnisotropic
author Martin Almquist <malmquist@stanford.edu>
date Mon, 29 Jun 2020 15:23:01 -0700
parents 15865fbda16e
children fe712a1fca3f
comparison
equal deleted inserted replaced
1281:01a0500de446 1291:a8e730db76e9
24 D % Total operator 24 D % Total operator
25 25
26 Dx, Dy % Physical derivatives 26 Dx, Dy % Physical derivatives
27 sigma % Cell matrix of physical stress operators 27 sigma % Cell matrix of physical stress operators
28 n_w, n_e, n_s, n_n % Physical normals 28 n_w, n_e, n_s, n_n % Physical normals
29 tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents
29 30
30 % Boundary operators in cell format, used for BC 31 % Boundary operators in cell format, used for BC
31 T_w, T_e, T_s, T_n 32 T_w, T_e, T_s, T_n
32 33
33 % Traction operators 34 % Traction operators
48 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary 49 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
49 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary 50 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
50 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary 51 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
51 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field 52 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
52 en_w, en_e, en_s, en_n % Act on vector field, return normal component 53 en_w, en_e, en_s, en_n % Act on vector field, return normal component
54 et_w, et_e, et_s, et_n % Act on vector field, return tangential component
53 55
54 % E{i}^T picks out component i 56 % E{i}^T picks out component i
55 E 57 E
56 58
57 % Elastic2dVariableAnisotropic object for reference domain 59 % Elastic2dVariableAnisotropic object for reference domain
278 obj.n_w = cell(2,1); 280 obj.n_w = cell(2,1);
279 obj.n_e = cell(2,1); 281 obj.n_e = cell(2,1);
280 obj.n_s = cell(2,1); 282 obj.n_s = cell(2,1);
281 obj.n_n = cell(2,1); 283 obj.n_n = cell(2,1);
282 284
285 % Compute normal and rotate (exactly!) 90 degrees counter-clockwise to get tangent
283 n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2))); 286 n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2)));
284 n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2))); 287 n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2)));
285 obj.n_w{1} = spdiag(n_w_1); 288 obj.n_w{1} = spdiag(n_w_1);
286 obj.n_w{2} = spdiag(n_w_2); 289 obj.n_w{2} = spdiag(n_w_2);
290 obj.tangent_w = {-obj.n_w{2}, obj.n_w{1}};
287 291
288 n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2))); 292 n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2)));
289 n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2))); 293 n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2)));
290 obj.n_e{1} = spdiag(n_e_1); 294 obj.n_e{1} = spdiag(n_e_1);
291 obj.n_e{2} = spdiag(n_e_2); 295 obj.n_e{2} = spdiag(n_e_2);
296 obj.tangent_e = {-obj.n_e{2}, obj.n_e{1}};
292 297
293 n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2))); 298 n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2)));
294 n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2))); 299 n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2)));
295 obj.n_s{1} = spdiag(n_s_1); 300 obj.n_s{1} = spdiag(n_s_1);
296 obj.n_s{2} = spdiag(n_s_2); 301 obj.n_s{2} = spdiag(n_s_2);
302 obj.tangent_s = {-obj.n_s{2}, obj.n_s{1}};
297 303
298 n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2))); 304 n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2)));
299 n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2))); 305 n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2)));
300 obj.n_n{1} = spdiag(n_n_1); 306 obj.n_n{1} = spdiag(n_n_1);
301 obj.n_n{2} = spdiag(n_n_2); 307 obj.n_n{2} = spdiag(n_n_2);
308 obj.tangent_n = {-obj.n_n{2}, obj.n_n{1}};
302 309
303 % Operators that extract the normal component 310 % Operators that extract the normal component
304 obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')'; 311 obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
305 obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')'; 312 obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')';
306 obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')'; 313 obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
307 obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')'; 314 obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
315
316 % Operators that extract the tangential component
317 obj.et_w = (obj.tangent_w{1}*obj.e1_w' + obj.tangent_w{2}*obj.e2_w')';
318 obj.et_e = (obj.tangent_e{1}*obj.e1_e' + obj.tangent_e{2}*obj.e2_e')';
319 obj.et_s = (obj.tangent_s{1}*obj.e1_s' + obj.tangent_s{2}*obj.e2_s')';
320 obj.et_n = (obj.tangent_n{1}*obj.e1_n' + obj.tangent_n{2}*obj.e2_n')';
308 321
309 % Stress operators 322 % Stress operators
310 sigma = cell(dim, dim); 323 sigma = cell(dim, dim);
311 D1 = {obj.Dx, obj.Dy}; 324 D1 = {obj.Dx, obj.Dy};
312 E = obj.E; 325 E = obj.E;
352 % In this way, we can specify one BC at a time even though the SATs depend on all BC. 365 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
353 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning) 366 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
354 default_arg('tuning', 1.0); 367 default_arg('tuning', 1.0);
355 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' ); 368 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
356 369
357 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning); 370 component = bc{1};
358
359 type = bc{2}; 371 type = bc{2};
372
373 switch component
374
375 % If conditions on Cartesian components
376 case {1,2}
377 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
378
379 % If conditions on normal or tangential components
380 case {'n', 't'}
381
382 switch component
383 case 'n'
384 en = obj.getBoundaryOperator('en', boundary);
385 case 't'
386 en = obj.getBoundaryOperator('et', boundary);
387 end
388 e1 = obj.getBoundaryOperator('e1', boundary);
389
390 bc = {1, type};
391 [c1, p1] = obj.refObj.boundary_condition(boundary, bc, tuning);
392 bc = {2, type};
393 c2 = obj.refObj.boundary_condition(boundary, bc, tuning);
394
395 switch type
396 case {'F','f','Free','free','traction','Traction','t','T'}
397 closure = en*en'*(c1+c2);
398 penalty = en*e1'*p1;
399 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
400 error('Not implemented');
401 end
402
403 end
360 404
361 switch type 405 switch type
362 case {'F','f','Free','free','traction','Traction','t','T'} 406 case {'F','f','Free','free','traction','Traction','t','T'}
407
363 s = obj.(['s_' boundary]); 408 s = obj.(['s_' boundary]);
364 s = spdiag(s); 409 s = spdiag(s);
365 penalty = penalty*s; 410 penalty = penalty*s;
411
366 end 412 end
367 end 413 end
368 414
369 % type Struct that specifies the interface coupling. 415 % type Struct that specifies the interface coupling.
370 % Fields: 416 % Fields:
402 448
403 % Returns the boundary operator op for the boundary specified by the string boundary. 449 % Returns the boundary operator op for the boundary specified by the string boundary.
404 % op -- string 450 % op -- string
405 function o = getBoundaryOperator(obj, op, boundary) 451 function o = getBoundaryOperator(obj, op, boundary)
406 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 452 assertIsMember(boundary, {'w', 'e', 's', 'n'})
407 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'}) 453 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et'})
408 454
409 o = obj.([op, '_', boundary]); 455 o = obj.([op, '_', boundary]);
410 456
411 end 457 end
412 458