annotate +scheme/Divergence.m @ 1207:05a01f77d0e3 feature/poroelastic

Swap indices and fix bug in ElasticVariableAnisotropic
author Martin Almquist <malmquist@stanford.edu>
date Fri, 20 Sep 2019 14:57:02 -0700
parents c2d281633e14
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1119
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
1 classdef Divergence < scheme.Scheme
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 % Approximates the divergence
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 % Interface and boundary condition methods are just dummies
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 properties
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 m % Number of points in each direction, possibly a vector
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 h % Grid spacing
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10 grid
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 dim
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13 order % Order of accuracy for the approximation
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 D
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16 D1
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17 H
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 methods
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 function obj = Divergence(g, order, opSet)
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25 dim = 2;
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 m = g.size();
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28 m_tot = g.N();
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 h = g.scaling();
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31 lim = g.lim;
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32 if isempty(lim)
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 x = g.x;
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 lim = cell(length(x),1);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35 for i = 1:length(x)
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36 lim{i} = {min(x{i}), max(x{i})};
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 % 1D operators
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 ops = cell(dim,1);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 for i = 1:dim
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 ops{i} = opSet{i}(m(i), lim{i}, order);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 I = cell(dim,1);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 D1 = cell(dim,1);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 for i = 1:dim
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 I{i} = speye(m(i));
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51 D1{i} = ops{i}.D1;
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54 %====== Assemble full operators ========
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56 % D1
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 obj.D1{1} = kron(D1{1},I{2});
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 obj.D1{2} = kron(I{1},D1{2});
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60 I_dim = speye(dim, dim);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62 % E{i}^T picks out component i.
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63 E = cell(dim,1);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64 I = speye(m_tot,m_tot);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
65 for i = 1:dim
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
66 e = sparse(dim,1);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
67 e(i) = 1;
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
68 E{i} = kron(I,e);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
69 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
70
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
71 Div = sparse(m_tot, dim*m_tot);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
72 for i = 1:dim
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
73 Div = Div + obj.D1{i}*E{i}';
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
74 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
75 obj.D = Div;
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
76 obj.H = [];
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
77 %=========================================%'
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
78
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
79 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
80
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
81
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
82 % Closure functions return the operators applied to the own domain to close the boundary
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
83 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
84 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
85 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
86 % on the first component. Can also be e.g.
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
87 % {'normal', 'd'} or {'tangential', 't'} for conditions on
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
88 % tangential/normal component.
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
89 % data is a function returning the data that should be applied at the boundary.
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90 % neighbour_scheme is an instance of Scheme that should be interfaced to.
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 % neighbour_boundary is a string specifying which boundary to interface to.
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92 function [closure, penalty] = boundary_condition(obj, boundary, bc)
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93 error('Not implemented')
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96 % type Struct that specifies the interface coupling.
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97 % Fields:
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
98 % -- tuning: penalty strength, defaults to 1.2
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
99 % -- interpolation: type of interpolation, default 'none'
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
100 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
1124
c2d281633e14 Divergence and Gradient do not work for multi-block. Not sure how to fix, left error message for now.
Martin Almquist <malmquist@stanford.edu>
parents: 1119
diff changeset
101 error('Not implemented')
1119
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
102 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
103
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
104 function N = size(obj)
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
105 N = obj.dim*prod(obj.m);
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
106 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
107 end
8984b12feba6 Add diffOp (currently scheme) for Divergence, not sure if this belongs in +scheme.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
108 end