Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_lonely_4_min_boundary_points.m @ 325:72468bc9b63f feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Renamed some operator implementations.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 26 Sep 2016 09:55:16 +0200
parents +sbp/+implementations/d4_variable_4_min_boundary_points.m@c0cbffcf6513
children b19e142fcae1
comparison
equal deleted inserted replaced
324:c0cbffcf6513 325:72468bc9b63f
1 function [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_4_min_boundary_points(m,h)
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %%% 4:de ordn. SBP Finita differens %%%
4 %%% operatorer framtagna av Mark Carpenter %%%
5 %%% %%%
6 %%% H (Normen) %%%
7 %%% D1=H^(-1)Q (approx f?rsta derivatan) %%%
8 %%% D2 (approx andra derivatan) %%%
9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10 %H?r med endast 4 randpunkter
11
12
13 BP = 4;
14 if(m<2*BP)
15 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
16 end
17
18
19 % Norm
20 Hv = ones(m,1);
21 Hv(1:4) = [17/48 59/48 43/48 49/48];
22 Hv(m-3:m) = rot90(Hv(1:4),2);
23 Hv = h*Hv;
24 H = spdiag(Hv, 0);
25 HI = spdiag(1./Hv, 0);
26
27
28 % Boundary operators
29 e_l = sparse(m,1);
30 e_l(1) = 1;
31 e_r = rot90(e_l, 2);
32
33 d1_l = sparse(m,1);
34 d1_l(1:4) = 1/h*[-11/6 3 -3/2 1/3];
35 d1_r = -rot90(d1_l);
36
37 d2_l = sparse(m,1);
38 d2_l(1:4) = 1/h^2*[2 -5 4 -1];
39 d2_r = rot90(d2_l, 2);
40
41 d3_l = sparse(m,1);
42 d3_l(1:4) = 1/h^3*[-1 3 -3 1];
43 d3_r = -rot90(d3_l, 2);
44
45
46 % First derivative
47 stencil = [1/12 -2/3 0 2/3 -1/12];
48 diags = [-1 0 1];
49
50 Q_U = [
51 0 0.59e2/0.96e2 -0.1e1/0.12e2 -0.1e1/0.32e2;
52 -0.59e2/0.96e2 0 0.59e2/0.96e2 0;
53 0.1e1/0.12e2 -0.59e2/0.96e2 0 0.59e2/0.96e2;
54 0.1e1/0.32e2 0 -0.59e2/0.96e2 0;
55 ];
56
57 Q = stripeMatrix(stencil, diags, m);
58 Q(1:4,1:4)=Q_U;
59 Q(m-3:m,m-3:m) = -rot90(Q_U, 2);
60
61 D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r');
62
63 % Fourth derivative
64 stencil = [-1/6, 2, -13/2, 28/3, -13/2, 2, -1/6];
65 diags = -3:3;
66 M4 = stripeMatrix(stencil, diags, m);
67
68 M4_U=[
69 0.8e1/0.3e1 -0.37e2/0.6e1 0.13e2/0.3e1 -0.5e1/0.6e1;
70 -0.37e2/0.6e1 0.47e2/0.3e1 -13 0.11e2/0.3e1;
71 0.13e2/0.3e1 -13 0.44e2/0.3e1 -0.47e2/0.6e1;
72 -0.5e1/0.6e1 0.11e2/0.3e1 -0.47e2/0.6e1 0.29e2/0.3e1;
73 ];
74
75
76 M4(1:4,1:4) = M4_U;
77 M4(m-3:m,m-3:m) = rot90(M4_U, 2);
78 M4 = 1/h^3*M4;
79
80 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r');
81 end