Mercurial > repos > public > sbplib

comparison +scheme/Gradient.m @ 1121:dd25184e71ec feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add diffOp for gradient
author Martin Almquist <malmquist@stanford.edu>
date Sat, 11 May 2019 15:11:04 -0700
parents
children
comparison
equal deleted inserted replaced
1120:7963a9cab637 1121:dd25184e71ec
1 classdef Gradient < scheme.Scheme
2
3 % Approximates the divergence
4 % Interface and boundary condition methods are just dummies
5
6 properties
7 m % Number of points in each direction, possibly a vector
8 h % Grid spacing
9
10 grid
11 dim
12
13 order % Order of accuracy for the approximation
14
15 D
16 D1
17 H
18 end
19
20 methods
21
22 function obj = Gradient(g, order, opSet)
23 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
24
25 dim = 2;
26
27 m = g.size();
28 m_tot = g.N();
29
30 h = g.scaling();
31 lim = g.lim;
32 if isempty(lim)
33 x = g.x;
34 lim = cell(length(x),1);
35 for i = 1:length(x)
36 lim{i} = {min(x{i}), max(x{i})};
37 end
38 end
39
40 % 1D operators
41 ops = cell(dim,1);
42 for i = 1:dim
43 ops{i} = opSet{i}(m(i), lim{i}, order);
44 end
45
46 I = cell(dim,1);
47 D1 = cell(dim,1);
48
49 for i = 1:dim
50 I{i} = speye(m(i));
51 D1{i} = ops{i}.D1;
52 end
53
54 %====== Assemble full operators ========
55
56 % D1
57 obj.D1{1} = kron(D1{1},I{2});
58 obj.D1{2} = kron(I{1},D1{2});
59
60 I_dim = speye(dim, dim);
61
62 % E{i}^T picks out component i.
63 E = cell(dim,1);
64 I = speye(m_tot,m_tot);
65 for i = 1:dim
66 e = sparse(dim,1);
67 e(i) = 1;
68 E{i} = kron(I,e);
69 end
70
71 Grad = sparse(dim*m_tot, m_tot);
72 for i = 1:dim
73 Grad = Grad + E{i}*obj.D1{i};
74 end
75 obj.D = Grad;
76 obj.H = [];
77 %=========================================%'
78
79 end
80
81
82 % Closure functions return the operators applied to the own domain to close the boundary
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.
84 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
85 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
86 % on the first component. Can also be e.g.
87 % {'normal', 'd'} or {'tangential', 't'} for conditions on
88 % tangential/normal component.
89 % data is a function returning the data that should be applied at the boundary.
90 % neighbour_scheme is an instance of Scheme that should be interfaced to.
91 % neighbour_boundary is a string specifying which boundary to interface to.
92 function [closure, penalty] = boundary_condition(obj, boundary, bc)
93 error('Not implemented')
94 end
95
96 % type Struct that specifies the interface coupling.
97 % Fields:
98 % -- tuning: penalty strength, defaults to 1.2
99 % -- interpolation: type of interpolation, default 'none'
100 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
101
102 [m, n] = size(obj.D);
103 closure = sparse(m, n);
104
105 [m, n] = size(neighbour_scheme.D);
106 penalty = sparse(m, n);
107 end
108
109 function N = size(obj)
110 N = obj.dim*prod(obj.m);
111 end
112 end
113 end