Mercurial > repos > public > sbplib

comparison +scheme/elasticVariable.m @ 685:b035902869a8 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make scheme/elasticVariable work with different opSets in different dim to facilitate periodic in some dim.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 08 Feb 2018 16:43:43 -0800
parents 247b58a4dbe8
children
comparison
equal deleted inserted replaced
684:fcf004066ea9 685:b035902869a8
1 classdef elasticVariable < scheme.Scheme 1 classdef elasticVariable < scheme.Scheme
2 2
3 % Discretizes the elastic wave equation: 3 % Discretizes the elastic wave equation:
4 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i 4 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i
5 % opSet should be cell array of opSets, one per dimension. This
6 % is useful if we have periodic BC in one direction.
5 7
6 properties 8 properties
7 m % Number of points in each direction, possibly a vector 9 m % Number of points in each direction, possibly a vector
8 h % Grid spacing 10 h % Grid spacing
9 11
44 end 46 end
45 47
46 methods 48 methods
47 49
48 function obj = elasticVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet) 50 function obj = elasticVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet)
49 default_arg('opSet',@sbp.D2Variable); 51 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
50 default_arg('lambda_fun', @(x,y) 0*x+1); 52 default_arg('lambda_fun', @(x,y) 0*x+1);
51 default_arg('mu_fun', @(x,y) 0*x+1); 53 default_arg('mu_fun', @(x,y) 0*x+1);
52 default_arg('rho_fun', @(x,y) 0*x+1); 54 default_arg('rho_fun', @(x,y) 0*x+1);
53 dim = 2; 55 dim = 2;
54 56
59 rho = grid.evalOn(g, rho_fun); 61 rho = grid.evalOn(g, rho_fun);
60 m = g.size(); 62 m = g.size();
61 m_tot = g.N(); 63 m_tot = g.N();
62 64
63 h = g.scaling(); 65 h = g.scaling();
64 L = (m-1).*h; 66 lim = g.lim;
65 67
66 % 1D operators 68 % 1D operators
67 ops = cell(dim,1); 69 ops = cell(dim,1);
68 for i = 1:dim 70 for i = 1:dim
69 ops{i} = opSet(m(i), {0, L(i)}, order); 71 ops{i} = opSet{i}(m(i), lim{i}, order);
70 end 72 end
71 73
72 % Borrowing constants 74 % Borrowing constants
73 beta = ops{1}.borrowing.R.delta_D; 75 for i = 1:dim
74 obj.H11 = ops{1}.borrowing.H11; 76 beta = ops{i}.borrowing.R.delta_D;
75 obj.phi = beta/obj.H11; 77 obj.H11{i} = ops{i}.borrowing.H11;
76 obj.gamma = ops{1}.borrowing.M.d1; 78 obj.phi{i} = beta/obj.H11{i};
79 obj.gamma{i} = ops{i}.borrowing.M.d1;
80 end
77 81
78 I = cell(dim,1); 82 I = cell(dim,1);
79 D1 = cell(dim,1); 83 D1 = cell(dim,1);
80 D2 = cell(dim,1); 84 D2 = cell(dim,1);
81 H = cell(dim,1); 85 H = cell(dim,1);
310 314
311 % Dirichlet boundary condition 315 % Dirichlet boundary condition
312 case {'D','d','dirichlet','Dirichlet'} 316 case {'D','d','dirichlet','Dirichlet'}
313 317
314 tuning = 1.2; 318 tuning = 1.2;
315 phi = obj.phi; 319 phi = obj.phi{j};
316 h = obj.h(j); 320 h = obj.h(j);
317 h11 = obj.H11*h; 321 h11 = obj.H11{j}*h;
318 gamma = obj.gamma; 322 gamma = obj.gamma{j};
319 323
320 a_lambda = dim/h11 + 1/(h11*phi); 324 a_lambda = dim/h11 + 1/(h11*phi);
321 a_mu_i = 2/(gamma*h); 325 a_mu_i = 2/(gamma*h);
322 a_mu_ij = 2/h11 + 1/(h11*phi); 326 a_mu_ij = 2/h11 + 1/(h11*phi);
323 327