Mercurial > repos > public > sbplib
annotate +scheme/elasticShearVariable.m @ 667:ed853945ee99 feature/poroelastic
Make free BC work with data. Bugfix domain size.
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Wed, 20 Dec 2017 06:52:17 +0100 |
parents | 8e6dfd22fc59 |
children | 17e62551cdc2 |
rev | line source |
---|---|
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
1 classdef elasticShearVariable < scheme.Scheme |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
2 properties |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
3 m % Number of points in each direction, possibly a vector |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
4 h % Grid spacing |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
5 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
6 grid |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
7 dim |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
8 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
9 order % Order accuracy for the approximation |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
10 |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
11 A % Variable coefficient lambda of the operator (as diagonal matrix here) |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
12 RHO % Density (as diagonal matrix here) |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
13 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
14 D % Total operator |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
15 D1 % First derivatives |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
16 D2 % Second derivatives |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
17 H, Hi % Inner products |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
18 e_l, e_r |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
19 d1_l, d1_r % Normal derivatives at the boundary |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
20 E % E{i}^T picks out component i |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
21 |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
22 H_boundary % Boundary inner products |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
23 |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
24 A_boundary_l % Variable coefficient at boundaries |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
25 A_boundary_r % |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
26 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
27 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
28 methods |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
29 % Implements the shear part of the elastic wave equation, i.e. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
30 % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
31 % where a = lambda. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
32 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
33 function obj = elasticShearVariable(g ,order, a_fun, rho_fun, opSet) |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
34 default_arg('opSet',@sbp.D2Variable); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
35 default_arg('a_fun', @(x,y) 0*x+1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
36 default_arg('rho_fun', @(x,y) 0*x+1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
37 dim = 2; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
38 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
39 assert(isa(g, 'grid.Cartesian')) |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
40 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
41 a = grid.evalOn(g, a_fun); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
42 rho = grid.evalOn(g, rho_fun); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
43 m = g.size(); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
44 m_tot = g.N(); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
45 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
46 h = g.scaling(); |
667
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
47 L = (m-1).*h; |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
48 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
49 % 1D operators |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
50 ops = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
51 for i = 1:dim |
667
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
52 ops{i} = opSet(m(i), {0, L(i)}, order); |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
53 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
54 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
55 I = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
56 D1 = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
57 D2 = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
58 H = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
59 Hi = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
60 e_l = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
61 e_r = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
62 d1_l = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
63 d1_r = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
64 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
65 for i = 1:dim |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
66 I{i} = speye(m(i)); |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
67 D1{i} = ops{i}.D1; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
68 D2{i} = ops{i}.D2; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
69 H{i} = ops{i}.H; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
70 Hi{i} = ops{i}.HI; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
71 e_l{i} = ops{i}.e_l; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
72 e_r{i} = ops{i}.e_r; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
73 d1_l{i} = ops{i}.d1_l; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
74 d1_r{i} = ops{i}.d1_r; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
75 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
76 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
77 %====== Assemble full operators ======== |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
78 A = spdiag(a); |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
79 obj.A = A; |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
80 RHO = spdiag(rho); |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
81 obj.RHO = RHO; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
82 |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
83 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
84 obj.D1 = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
85 obj.D2 = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
86 obj.e_l = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
87 obj.e_r = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
88 obj.d1_l = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
89 obj.d1_r = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
90 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
91 % D1 |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
92 obj.D1{1} = kron(D1{1},I{2}); |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
93 obj.D1{2} = kron(I{1},D1{2}); |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
94 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
95 % Boundary operators |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
96 obj.e_l{1} = kron(e_l{1},I{2}); |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
97 obj.e_l{2} = kron(I{1},e_l{2}); |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
98 obj.e_r{1} = kron(e_r{1},I{2}); |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
99 obj.e_r{2} = kron(I{1},e_r{2}); |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
100 |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
101 obj.d1_l{1} = kron(d1_l{1},I{2}); |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
102 obj.d1_l{2} = kron(I{1},d1_l{2}); |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
103 obj.d1_r{1} = kron(d1_r{1},I{2}); |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
104 obj.d1_r{2} = kron(I{1},d1_r{2}); |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
105 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
106 % D2 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
107 for i = 1:dim |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
108 obj.D2{i} = sparse(m_tot); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
109 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
110 ind = grid.funcToMatrix(g, 1:m_tot); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
111 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
112 for i = 1:m(2) |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
113 D = D2{1}(a(ind(:,i))); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
114 p = ind(:,i); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
115 obj.D2{1}(p,p) = D; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
116 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
117 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
118 for i = 1:m(1) |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
119 D = D2{2}(a(ind(i,:))); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
120 p = ind(i,:); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
121 obj.D2{2}(p,p) = D; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
122 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
123 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
124 % Quadratures |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
125 obj.H = kron(H{1},H{2}); |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
126 obj.Hi = inv(obj.H); |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
127 obj.H_boundary = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
128 obj.H_boundary{1} = H{2}; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
129 obj.H_boundary{2} = H{1}; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
130 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
131 % Boundary coefficient matrices and quadratures |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
132 obj.A_boundary_l = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
133 obj.A_boundary_r = cell(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
134 for i = 1:dim |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
135 obj.A_boundary_l{i} = obj.e_l{i}'*A*obj.e_l{i}; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
136 obj.A_boundary_r{i} = obj.e_r{i}'*A*obj.e_r{i}; |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
137 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
138 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
139 % E{i}^T picks out component i. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
140 E = cell(dim,1); |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
141 I = speye(m_tot,m_tot); |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
142 for i = 1:dim |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
143 e = sparse(dim,1); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
144 e(i) = 1; |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
145 E{i} = kron(I,e); |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
146 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
147 obj.E = E; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
148 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
149 % Differentiation matrix D (without SAT) |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
150 D2 = obj.D2; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
151 D1 = obj.D1; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
152 D = sparse(dim*m_tot,dim*m_tot); |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
153 d = @kroneckerDelta; % Kronecker delta |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
154 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
155 for i = 1:dim |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
156 for j = 1:dim |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
157 D = D + E{i}*inv(RHO)*( d(i,j)*D2{i}*E{j}' +... |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
158 db(i,j)*D1{j}*A*D1{i}*E{j}' + ... |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
159 D2{j}*E{i}' ... |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
160 ); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
161 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
162 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
163 obj.D = D; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
164 %=========================================% |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
165 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
166 % Misc. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
167 obj.m = m; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
168 obj.h = h; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
169 obj.order = order; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
170 obj.grid = g; |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
171 obj.dim = dim; |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
172 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
173 % obj.gamm_u = h_u*ops_u.borrowing.M.d1; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
174 % obj.gamm_v = h_v*ops_v.borrowing.M.d1; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
175 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
176 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
177 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
178 % Closure functions return the operators applied to the own domain to close the boundary |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
179 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
180 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
181 % type is a string specifying the type of boundary condition if there are several. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
182 % data is a function returning the data that should be applied at the boundary. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
183 % neighbour_scheme is an instance of Scheme that should be interfaced to. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
184 % neighbour_boundary is a string specifying which boundary to interface to. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
185 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
186 default_arg('type','free'); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
187 default_arg('parameter', []); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
188 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
189 delta = @kroneckerDelta; % Kronecker delta |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
190 delta_b = @(i,j) 1-delta(i,j); % Logical not of Kronecker delta |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
191 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
192 % j is the coordinate direction of the boundary |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
193 % nj: outward unit normal component. |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
194 % nj = -1 for west, south, bottom boundaries |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
195 % nj = 1 for east, north, top boundaries |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
196 [j, nj] = obj.get_boundary_number(boundary); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
197 switch nj |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
198 case 1 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
199 e = obj.e_r; |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
200 d = obj.d1_r; |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
201 case -1 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
202 e = obj.e_l; |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
203 d = obj.d1_l; |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
204 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
205 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
206 E = obj.E; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
207 Hi = obj.Hi; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
208 H_gamma = obj.H_boundary{j}; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
209 A = obj.A; |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
210 RHO = obj.RHO; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
211 D1 = obj.D1; |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
212 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
213 switch type |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
214 % Dirichlet boundary condition |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
215 case {'D','d','dirichlet'} |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
216 error('Dirichlet not implemented') |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
217 tuning = 1.2; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
218 % tuning = 20.2; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
219 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
220 b1 = gamm*obj.lambda./obj.a11.^2; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
221 b2 = gamm*obj.lambda./obj.a22.^2; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
222 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
223 tau1 = tuning * spdiag(-1./b1 - 1./b2); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
224 tau2 = 1; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
225 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
226 tau = (tau1*e + tau2*d)*H_b; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
227 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
228 closure = obj.a*obj.Hi*tau*e'; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
229 penalty = -obj.a*obj.Hi*tau; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
230 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
231 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
232 % Free boundary condition |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
233 case {'F','f','Free','free'} |
667
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
234 closures = cell(obj.dim,1); |
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
235 penalties = cell(obj.dim,1); |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
236 % Loop over components |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
237 for i = 1:obj.dim |
667
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
238 closures{i} = E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma*(... |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
239 e{j}'*A*e{j}*d{j}'*E{i}' + ... |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
240 delta(i,j)*e{j}'*A*e{i}*d{i}'*E{j}' + ... |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
241 delta_b(i,j)*e{j}'*A*D1{i}*E{j}' ... |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
242 ); |
667
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
243 penalties{i} = - E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma; |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
244 end |
667
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
245 [rows, cols] = size(closures{1}); |
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
246 closure = sparse(rows, cols); |
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
247 for i = 1:obj.dim |
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
248 closure = closure + closures{i}; |
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
249 end |
ed853945ee99
Make free BC work with data. Bugfix domain size.
Martin Almquist <malmquist@stanford.edu>
parents:
664
diff
changeset
|
250 penalty = penalties; |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
251 |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
252 |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
253 % Unknown boundary condition |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
254 otherwise |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
255 error('No such boundary condition: type = %s',type); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
256 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
257 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
258 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
259 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
260 % u denotes the solution in the own domain |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
261 % v denotes the solution in the neighbour domain |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
262 tuning = 1.2; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
263 % tuning = 20.2; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
264 error('Interface not implemented'); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
265 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
266 |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
267 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
268 function [j, nj] = get_boundary_number(obj, boundary) |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
269 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
270 switch boundary |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
271 case {'w','W','west','West', 'e', 'E', 'east', 'East'} |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
272 j = 1; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
273 case {'s','S','south','South', 'n', 'N', 'north', 'North'} |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
274 j = 2; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
275 otherwise |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
276 error('No such boundary: boundary = %s',boundary); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
277 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
278 |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
279 switch boundary |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
280 case {'w','W','west','West','s','S','south','South'} |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
281 nj = -1; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
282 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
283 nj = 1; |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
284 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
285 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
286 |
664
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
287 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
288 function [return_op] = get_boundary_operator(obj, op, boundary) |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
289 |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
290 switch boundary |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
291 case {'w','W','west','West', 'e', 'E', 'east', 'East'} |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
292 j = 1; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
293 case {'s','S','south','South', 'n', 'N', 'north', 'North'} |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
294 j = 2; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
295 otherwise |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
296 error('No such boundary: boundary = %s',boundary); |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
297 end |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
298 |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
299 switch op |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
300 case 'e' |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
301 switch boundary |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
302 case {'w','W','west','West','s','S','south','South'} |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
303 return_op = obj.e_l{j}; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
304 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
305 return_op = obj.e_r{j}; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
306 end |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
307 case 'd' |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
308 switch boundary |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
309 case {'w','W','west','West','s','S','south','South'} |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
310 return_op = obj.d_l{j}; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
311 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
312 return_op = obj.d_r{j}; |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
313 end |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
314 otherwise |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
315 error(['No such operator: operatr = ' op]); |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
316 end |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
317 |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
318 end |
8e6dfd22fc59
Free BC now yields symmetric negative semidefinite discretization.
Martin Almquist <malmquist@stanford.edu>
parents:
663
diff
changeset
|
319 |
663
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
320 function N = size(obj) |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
321 N = prod(obj.m); |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
322 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
323 end |
b45ec2b28cc2
First implementation of elastic shear operator with free boundary BC.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
324 end |