1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
| /**
#A gnu-Octave / matlab code
To generate your own refinement attributes.
*/
% ax+by+cxy
disp 2d
clear all
A=4*[1 0;1 -1;0 -1; -1 -1;-1 0;-1 1;0 1; 1 1];
x=A(:,1);
y=A(:,2);
X(:,1)=x;
X(:,2)=y;
X(:,3)=(x.*y);
transpose(A)
C=inv(transpose(X)*X)*transpose(X);
C*192
%ax+by+cz+dxy+eyz+fxz+gxyz
disp 3dc
clear all
Locs=[ 1 0 0;1 -1 0;0 -1 0; -1 -1 0;-1 0 0;-1 1 0;0 1 0;1 1 0;
0 0 -1; 1 0 -1;1 -1 -1;0 -1 -1; -1 -1 -1;-1 0 -1;-1 1 -1;0 1 -1;1 1 -1;
0 0 1; 1 0 1;1 -1 1;0 -1 1; -1 -1 1;-1 0 1;-1 1 1;0 1 1;1 1 1];
A=4*Locs;
x=A(:,1);
y=A(:,2);
z=A(:,3);
X(:,1)=x;
X(:,2)=y;
X(:,3)=z;
X(:,4)=(x.*y);
X(:,5)=(z.*y);
X(:,6)=(x.*z);
X(:,7)=(x.*y.*z);
transpose(A)
C=inv(transpose(X)*X)*transpose(X);
C=C*8*9*8*8;
[a,~]=size(C);
for j=1:a
ind=find(C(j,:)~=0);
fprintf([num2str(round(abs(C(j,ind(1))))) '*']);
if (length(unique(abs(C(j,ind))))==1)
for g=1:length(ind);
if (C(j,ind(g))<0)
fprintf('-');
end
if (C(j,ind(g))>0 && g~=1)
fprintf('+');
end
fprintf(['s[' num2str(Locs(ind(g),1)) ',' num2str(Locs(ind(g),2)) ',' num2str(Locs(ind(g),3)) ']']);
end
clear ind
fprintf(';\n');
else
disp('Error! non unique weights');
end
end
%ax+by+cz+dxy+eyz+fxz+gxyz+hxx+iyy+jzz+kyxx+lzxx+mxyy+nzyy+oxzz+pzzy+qxxyy+rxxzz+syyzz+txxyyzz
%20 weights!
disp ('3d third order accurate')
clear all
Locs=[0 0 0; 1 0 0;1 -1 0;0 -1 0; -1 -1 0;-1 0 0;-1 1 0;0 1 0;1 1 0;
0 0 -1; 1 0 -1;1 -1 -1;0 -1 -1; -1 -1 -1;-1 0 -1;-1 1 -1;0 1 -1;1 1 -1;
0 0 1; 1 0 1;1 -1 1;0 -1 1; -1 -1 1;-1 0 1;-1 1 1;0 1 1;1 1 1;
2 0 0; -2 0 0; 0 2 0; 0 -2 0; 0 0 2; 0 0 -2];
A=4*Locs;
x=A(:,1);
y=A(:,2);
z=A(:,3);
% The terms
X(:,1)=x; %a
X(:,2)=y; %b
X(:,3)=z; %|
X(:,4)=(x.*y);
X(:,5)=(x.*z);
X(:,6)=(y.*z);
X(:,7)=(x.*y.*z);
X(:,8)=x.*x;
X(:,9)=y.*y;
X(:,10)=z.*z;
X(:,11)=X(:,8).*y;
X(:,12)=X(:,8).*z;
X(:,13)=X(:,9).*x;
X(:,14)=X(:,9).*z;
X(:,15)=X(:,10).*x;
X(:,16)=X(:,10).*y;
X(:,17)=X(:,8).*X(:,9);
X(:,18)=X(:,8).*X(:,10);
X(:,19)=X(:,9).*X(:,10); %|
X(:,20)=X(:,8).*X(:,9).*X(:,10); %t
X(:,21)=ones(size(x));
S{:}=['a';'b';'c';'d';'e';'f';'g';'h';'i';'j';'k';'l';'m';'n';'o';'p';'q';'r';'s';'t';'u'];
C=inv(transpose(X)*X)*transpose(X);
[a,~]=size(C);
C=round(C*100000,9);
for j=1:a
fprintf('double %c = ',S{1}(j));
ind=find(C(j,:)~=0);
wei=(unique(abs(C(j,ind))));
for h=1:length(wei);
ind = find(abs(C(j,:))==wei(h));
fprintf('%.15g*(',(wei(h)));
for g=1:length(ind);
if (C(j,ind(g))<0)
fprintf('-');
end
if (C(j,ind(g))>0 && g~=1)
fprintf('+');
end
fprintf(['s[' num2str(Locs(ind(g),1)) ',' num2str(Locs(ind(g),2)) ',' num2str(Locs(ind(g),3)) ']']);
end
clear ind
if h==length(wei)
fprintf(');\n');
else
fprintf(')+');
end
end
end
|
