Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropic.m @ 1213:43f1cd11e8e8 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add physical normals to AnisotropicCurvilinear
author Martin Almquist <malmquist@stanford.edu>
date Mon, 14 Oct 2019 13:54:50 -0700
parents 3d7faa2ca312
children e1d4cb8b5309
comparison
equal deleted inserted replaced
1212:5c5815af4b7a 1213:43f1cd11e8e8
22 C % Elastic stiffness tensor 22 C % Elastic stiffness tensor
23 23
24 D % Total operator 24 D % Total operator
25 25
26 Dx, Dy % Physical derivatives 26 Dx, Dy % Physical derivatives
27 n_w, n_e, n_s, n_n % Physical normals
27 28
28 % Boundary operators in cell format, used for BC 29 % Boundary operators in cell format, used for BC
29 T_w, T_e, T_s, T_n 30 T_w, T_e, T_s, T_n
30 31
31 % Traction operators 32 % Traction operators
45 % Boundary restriction operators 46 % Boundary restriction operators
46 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary 47 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
47 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary 48 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
48 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary 49 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
49 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field 50 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
51 en_w, en_e, en_s, en_n % Act on vector field, return normal component
50 52
51 % E{i}^T picks out component i 53 % E{i}^T picks out component i
52 E 54 E
53 55
54 % Elastic2dVariableAnisotropic object for reference domain 56 % Elastic2dVariableAnisotropic object for reference domain
229 obj.e_scalar_w = refObj.e_scalar_w; 231 obj.e_scalar_w = refObj.e_scalar_w;
230 obj.e_scalar_e = refObj.e_scalar_e; 232 obj.e_scalar_e = refObj.e_scalar_e;
231 obj.e_scalar_s = refObj.e_scalar_s; 233 obj.e_scalar_s = refObj.e_scalar_s;
232 obj.e_scalar_n = refObj.e_scalar_n; 234 obj.e_scalar_n = refObj.e_scalar_n;
233 235
234 e1_w = (obj.e_scalar_w'*obj.E{1}')'; 236 e1_w = obj.e1_w;
235 e1_e = (obj.e_scalar_e'*obj.E{1}')'; 237 e1_e = obj.e1_e;
236 e1_s = (obj.e_scalar_s'*obj.E{1}')'; 238 e1_s = obj.e1_s;
237 e1_n = (obj.e_scalar_n'*obj.E{1}')'; 239 e1_n = obj.e1_n;
238 240
239 e2_w = (obj.e_scalar_w'*obj.E{2}')'; 241 e2_w = obj.e2_w;
240 e2_e = (obj.e_scalar_e'*obj.E{2}')'; 242 e2_e = obj.e2_e;
241 e2_s = (obj.e_scalar_s'*obj.E{2}')'; 243 e2_s = obj.e2_s;
242 e2_n = (obj.e_scalar_n'*obj.E{2}')'; 244 e2_n = obj.e2_n;
243 245
244 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')'; 246 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
245 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')'; 247 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
246 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')'; 248 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
247 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')'; 249 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
253 255
254 obj.tau_w = (refObj.e_w'*obj.e1_w*obj.tau1_w')' + (refObj.e_w'*obj.e2_w*obj.tau2_w')'; 256 obj.tau_w = (refObj.e_w'*obj.e1_w*obj.tau1_w')' + (refObj.e_w'*obj.e2_w*obj.tau2_w')';
255 obj.tau_e = (refObj.e_e'*obj.e1_e*obj.tau1_e')' + (refObj.e_e'*obj.e2_e*obj.tau2_e')'; 257 obj.tau_e = (refObj.e_e'*obj.e1_e*obj.tau1_e')' + (refObj.e_e'*obj.e2_e*obj.tau2_e')';
256 obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')'; 258 obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')';
257 obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')'; 259 obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')';
260
261 % Physical normals
262 e_w = obj.e_scalar_w;
263 e_e = obj.e_scalar_e;
264 e_s = obj.e_scalar_s;
265 e_n = obj.e_scalar_n;
266
267 e_w_vec = obj.e_w;
268 e_e_vec = obj.e_e;
269 e_s_vec = obj.e_s;
270 e_n_vec = obj.e_n;
271
272 nu_w = [-1,0];
273 nu_e = [1,0];
274 nu_s = [0,-1];
275 nu_n = [0,1];
276
277 obj.n_w = cell(2,1);
278 obj.n_e = cell(2,1);
279 obj.n_s = cell(2,1);
280 obj.n_n = cell(2,1);
281
282 n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2)));
283 n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2)));
284 obj.n_w{1} = spdiag(n_w_1);
285 obj.n_w{2} = spdiag(n_w_2);
286
287 n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2)));
288 n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2)));
289 obj.n_e{1} = spdiag(n_e_1);
290 obj.n_e{2} = spdiag(n_e_2);
291
292 n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2)));
293 n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2)));
294 obj.n_s{1} = spdiag(n_s_1);
295 obj.n_s{2} = spdiag(n_s_2);
296
297 n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2)));
298 n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2)));
299 obj.n_n{1} = spdiag(n_n_1);
300 obj.n_n{2} = spdiag(n_n_2);
301
302 % Operators that extract the normal component
303 obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
304 obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')';
305 obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
306 obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
258 307
259 % Misc. 308 % Misc.
260 obj.refObj = refObj; 309 obj.refObj = refObj;
261 obj.m = refObj.m; 310 obj.m = refObj.m;
262 obj.h = refObj.h; 311 obj.h = refObj.h;
310 default_struct('type', defaultType); 359 default_struct('type', defaultType);
311 360
312 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type); 361 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
313 end 362 end
314 363
315 % Returns the component number that is the tangential/normal component
316 % at the specified boundary
317 function comp = getComponent(obj, comp_str, boundary)
318 assertIsMember(comp_str, {'normal', 'tangential'});
319 assertIsMember(boundary, {'w', 'e', 's', 'n'});
320
321 switch boundary
322 case {'w', 'e'}
323 switch comp_str
324 case 'normal'
325 comp = 1;
326 case 'tangential'
327 comp = 2;
328 end
329 case {'s', 'n'}
330 switch comp_str
331 case 'normal'
332 comp = 2;
333 case 'tangential'
334 comp = 1;
335 end
336 end
337 end
338
339 % Returns h11 for the boundary specified by the string boundary. 364 % Returns h11 for the boundary specified by the string boundary.
340 % op -- string 365 % op -- string
341 function h11 = getBorrowing(obj, boundary) 366 function h11 = getBorrowing(obj, boundary)
342 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 367 assertIsMember(boundary, {'w', 'e', 's', 'n'})
343 368
348 h11 = obj.refObj.h11{2}; 373 h11 = obj.refObj.h11{2};
349 end 374 end
350 end 375 end
351 376
352 % Returns the outward unit normal vector for the boundary specified by the string boundary. 377 % Returns the outward unit normal vector for the boundary specified by the string boundary.
353 function nu = getNormal(obj, boundary) 378 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
354 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 379 function n = getNormal(obj, boundary)
355 380 assertIsMember(boundary, {'w', 'e', 's', 'n'})
356 switch boundary 381
357 case 'w' 382 n = obj.(['n_' boundary]);
358 nu = [-1,0];
359 case 'e'
360 nu = [1,0];
361 case 's'
362 nu = [0,-1];
363 case 'n'
364 nu = [0,1];
365 end
366 end 383 end
367 384
368 % Returns the boundary operator op for the boundary specified by the string boundary. 385 % Returns the boundary operator op for the boundary specified by the string boundary.
369 % op -- string 386 % op -- string
370 function o = getBoundaryOperator(obj, op, boundary) 387 function o = getBoundaryOperator(obj, op, boundary)
371 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 388 assertIsMember(boundary, {'w', 'e', 's', 'n'})
372 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}) 389 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'})
373 390
374 switch op 391 o = obj.([op, '_', boundary]);
375 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
376 o = obj.([op, '_', boundary]);
377 end
378 392
379 end 393 end
380 394
381 % Returns the boundary operator op for the boundary specified by the string boundary. 395 % Returns the boundary operator op for the boundary specified by the string boundary.
382 % op -- string 396 % op -- string