1 /* Copyright (C) 2003-2007 CAMP
2 * Copyright (C) 2005 CSC - IT Center for Science Ltd.
3 * Please see the accompanying LICENSE file for further information. */
4
5 #include <stdlib.h>
6 #include "bmgs.h"
7
8
9 // Expansion coefficients for finite difference Laplacian. The numbers are
10 // from J. R. Chelikowsky et al., Phys. Rev. B 50, 11355 (1994):
11
bmgs_stencil(int ncoefs,const double * coefs,const long * offsets,int r,const long n[3])12 bmgsstencil bmgs_stencil(int ncoefs, const double* coefs, const long* offsets,
13 int r, const long n[3])
14 {
15 bmgsstencil stencil =
16 {ncoefs,
17 (double*)malloc(ncoefs * sizeof(double)),
18 (long*)malloc(ncoefs * sizeof(long)),
19 {n[0], n[1], n[2]},
20 {2 * r * (n[2] + 2 * r) * (n[1] + 2 * r),
21 2 * r * (n[2] + 2 * r),
22 2 * r}};
23 assert((stencil.coefs != NULL) && (stencil.offsets != NULL));
24 memcpy(stencil.coefs, coefs, ncoefs * sizeof(double));
25 memcpy(stencil.offsets, offsets, ncoefs * sizeof(long));
26 return stencil;
27 }
28
29
30 static const double laplace[4][5] =
31 {{-2.0, 1.0, 0.0, 0.0, 0.0},
32 {-5.0/2.0, 4.0/3.0, -1.0/12.0, 0.0, 0.0},
33 {-49.0/18.0, 3.0/2.0, -3.0/20.0, 1.0/90.0, 0.0},
34 {-205.0/72.0, 8.0/5.0, -1.0/5.0, 8.0/315.0, -1.0/560.0}};
35
36
bmgs_laplace(int k,double scale,const double h[3],const long n[3])37 bmgsstencil bmgs_laplace(int k, double scale,
38 const double h[3],
39 const long n[3])
40 {
41 int ncoefs = 3 * k - 2;
42 double* coefs = (double*)malloc(ncoefs * sizeof(double));
43 long* offsets = (long*)malloc(ncoefs * sizeof(long));
44 assert((coefs != NULL) && (offsets != NULL));
45 double f1 = 1.0 / (h[0] * h[0]);
46 double f2 = 1.0 / (h[1] * h[1]);
47 double f3 = 1.0 / (h[2] * h[2]);
48 int r = (k - 1) / 2; // range
49 double s[3] = {(n[2] + 2 * r) * (n[1] + 2 * r), n[2] + 2 * r, 1};
50 int m = 0;
51 for (int j = 1; j <= r; j++)
52 {
53 double c = scale * laplace[r - 1][j];
54 coefs[m] = c * f1; offsets[m++] = -j * s[0];
55 coefs[m] = c * f1; offsets[m++] = +j * s[0];
56 coefs[m] = c * f2; offsets[m++] = -j * s[1];
57 coefs[m] = c * f2; offsets[m++] = +j * s[1];
58 coefs[m] = c * f3; offsets[m++] = -j;
59 coefs[m] = c * f3; offsets[m++] = +j;
60 }
61 double c = scale * laplace[r - 1][0];
62 coefs[m] = c * (f1 + f2 + f3); offsets[m] = 0;
63 bmgsstencil stencil =
64 {ncoefs, coefs, offsets,
65 {n[0], n[1], n[2]},
66 {2 * r * (n[2] + 2 * r) * (n[1] + 2 * r),
67 2 * r * (n[2] + 2 * r),
68 2 * r}};
69 return stencil;
70 }
71
72
bmgs_mslaplaceA(double scale,const double h[3],const long n[3])73 bmgsstencil bmgs_mslaplaceA(double scale,
74 const double h[3],
75 const long n[3])
76 {
77 int ncoefs = 19;
78 double* coefs = (double*)malloc(ncoefs * sizeof(double));
79 long* offsets = (long*)malloc(ncoefs * sizeof(long));
80 assert((coefs != NULL) && (offsets != NULL));
81 double e[3] = {-scale / (12.0 * h[0] * h[0]),
82 -scale / (12.0 * h[1] * h[1]),
83 -scale / (12.0 * h[2] * h[2])};
84 double f = -16.0 * (e[0] + e[1] + e[2]);
85 double g[3] = {10.0 * e[0] + 0.125 * f,
86 10.0 * e[1] + 0.125 * f,
87 10.0 * e[2] + 0.125 * f};
88 double s[3] = {(n[2] + 2) * (n[1] + 2), n[2] + 2, 1};
89 int m = 0;
90 coefs[m] = f;
91 offsets[m++] = 0;
92 for (int j = -1; j <= 1; j += 2)
93 {
94 coefs[m] = g[0];
95 offsets[m++] = j * s[0];
96 coefs[m] = g[1];
97 offsets[m++] = j * s[1];
98 coefs[m] = g[2];
99 offsets[m++] = j * s[2];
100 }
101 for (int j = -1; j <= 1; j += 2)
102 for (int k = -1; j <= 1; j += 2)
103 {
104 coefs[m] = e[1] + e[2];
105 offsets[m++] = -j * s[1] - k * s[2];
106 coefs[m] = e[0] + e[2];
107 offsets[m++] = -j * s[0] - k * s[2];
108 coefs[m] = e[0] + e[1];
109 offsets[m++] = -j * s[0] - k * s[1];
110 }
111 bmgsstencil stencil =
112 {ncoefs, coefs, offsets,
113 {n[0], n[1], n[2]},
114 {2 * s[0], 2 * s[1], 2}};
115 return stencil;
116 }
117
118
bmgs_mslaplaceB(const long n[3])119 bmgsstencil bmgs_mslaplaceB(const long n[3])
120 {
121 int ncoefs = 7;
122 double* coefs = (double*)malloc(ncoefs * sizeof(double));
123 long* offsets = (long*)malloc(ncoefs * sizeof(long));
124 assert((coefs != NULL) && (offsets != NULL));
125 double s[3] = {(n[2] + 2) * (n[1] + 2), n[2] + 2, 1};
126 int k = 0;
127 coefs[k] = 0.5;
128 offsets[k++] = 0;
129 for (int j = -1; j <= 1; j += 2)
130 {
131 coefs[k] = 1.0 / 12.0;
132 offsets[k++] = j * s[0];
133 coefs[k] = 1.0 / 12.0;
134 offsets[k++] = j * s[1];
135 coefs[k] = 1.0 / 12.0;
136 offsets[k++] = j * s[2];
137 }
138 bmgsstencil stencil =
139 {ncoefs, coefs, offsets,
140 {n[0], n[1], n[2]},
141 {2 * s[0], 2 * s[1], 2}};
142 return stencil;
143 }
144
145
bmgs_gradient(int k,int i,double h,const long n[3])146 bmgsstencil bmgs_gradient(int k, int i, double h,
147 const long n[3])
148 {
149 int ncoefs = k - 1;
150 double* coefs = (double*)malloc(ncoefs * sizeof(double));
151 long* offsets = (long*)malloc(ncoefs * sizeof(long));
152 assert((coefs != NULL) && (offsets != NULL));
153 int r = 1;
154 double s[3] = {(n[2] + 2 * r) * (n[1] + 2 * r), n[2] + 2 * r, 1};
155 double c = 0.5 / h;
156 coefs[0] = +c; offsets[0] = +s[i];
157 coefs[1] = -c; offsets[1] = -s[i];
158 bmgsstencil stencil =
159 {ncoefs, coefs, offsets,
160 {n[0], n[1], n[2]},
161 {2 * r * s[0], 2 * r * s[1], 2 * r}};
162 return stencil;
163 }
164
bmgs_deletestencil(bmgsstencil * stencil)165 void bmgs_deletestencil(bmgsstencil* stencil)
166 {
167 free(stencil->coefs);
168 free(stencil->offsets);
169 }
170