Peano 4
Loading...
Searching...
No Matches
TP_Utilities.cpp
Go to the documentation of this file.
1/* TwoPunctures: File "utilities.c"*/
2/* these functions have no dependency on the TP parameters */
3
4#include <math.h>
5#include <stdio.h>
6#include <stddef.h>
7#include <stdlib.h>
8#include <assert.h>
10
11namespace TP {
12namespace Utilities {
13
14/*---------------------------------------------------------------------------*/
15int *
16ivector (long nl, long nh)
17/* allocate an int vector with subscript range v[nl..nh] */
18{
19 int *retval;
20
21 retval = new int [(nh-nl+1)];
22 if(retval == NULL)
23 throw error ("allocation failure in ivector()");
24
25 return retval - nl;
26}
27
28/*---------------------------------------------------------------------------*/
29double *
30dvector (long nl, long nh)
31/* allocate a double vector with subscript range v[nl..nh] */
32{
33 double *retval;
34
35 retval = new double [(nh-nl+1)];
36 if(retval == NULL)
37 throw error ("allocation failure in dvector()");
38
39 return retval - nl;
40}
41
42/*---------------------------------------------------------------------------*/
43int **
44imatrix (long nrl, long nrh, long ncl, long nch)
45/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
46{
47 int **retval;
48
49 retval = new int * [(nrh-nrl+1)];
50 if(retval == NULL)
51 throw error ("allocation failure (1) in imatrix()");
52
53 /* get all memory for the matrix in on chunk */
54 retval[0] = new int [(nrh-nrl+1)*(nch-ncl+1)];
55 if(retval[0] == NULL)
56 throw error ("allocation failure (2) in imatrix()");
57
58 /* apply column and row offsets */
59 retval[0] -= ncl;
60 retval -= nrl;
61
62 /* slice chunk into rows */
63 long width = (nch-ncl+1);
64 for(long i = nrl+1 ; i <= nrh ; i++)
65 retval[i] = retval[i-1] + width;
66 assert(retval[nrh]-retval[nrl] == (nrh-nrl)*width);
67
68 return retval;
69}
70
71/*---------------------------------------------------------------------------*/
72double **
73dmatrix (long nrl, long nrh, long ncl, long nch)
74/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
75{
76 double **retval;
77
78 retval = new double * [(nrh-nrl+1)];
79 if(retval == NULL)
80 throw error ("allocation failure (1) in dmatrix()");
81
82 /* get all memory for the matrix in on chunk */
83 retval[0] = new double [(nrh-nrl+1)*(nch-ncl+1)];
84 if(retval[0] == NULL)
85 throw error ("allocation failure (2) in dmatrix()");
86
87 /* apply column and row offsets */
88 retval[0] -= ncl;
89 retval -= nrl;
90
91 /* slice chunk into rows */
92 long width = (nch-ncl+1);
93 for(long i = nrl+1 ; i <= nrh ; i++)
94 retval[i] = retval[i-1] + width;
95 assert(retval[nrh]-retval[nrl] == (nrh-nrl)*width);
96
97 return retval;
98}
99
100/*---------------------------------------------------------------------------*/
101double ***
102d3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
103/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
104{
105 double ***retval;
106
107 /* get memory for index structures */
108 retval = new double ** [(nrh-nrl+1)];
109 if(retval == NULL)
110 throw error ("allocation failure (1) in d3tensor()");
111
112 retval[0] = new double * [(nrh-nrl+1)*(nch-ncl+1)];
113 if(retval[0] == NULL)
114 throw error ("allocation failure (2) in d3tensor()");
115
116 /* get all memory for the tensor in on chunk */
117 retval[0][0] = new double [(nrh-nrl+1)*(nch-ncl+1)*(ndh-ndl+1)];
118 if(retval[0][0] == NULL)
119 throw error ("allocation failure (3) in d3tensor()");
120
121 /* apply all offsets */
122 retval[0][0] -= ndl;
123 retval[0] -= ncl;
124 retval -= nrl;
125
126 /* slice chunk into rows and columns */
127 long width = (nch-ncl+1);
128 long depth = (ndh-ndl+1);
129 for(long j = ncl+1 ; j <= nch ; j++) { /* first row of columns */
130 retval[nrl][j] = retval[nrl][j-1] + depth;
131 }
132 assert(retval[nrl][nch]-retval[nrl][ncl] == (nch-ncl)*depth);
133 for(long i = nrl+1 ; i <= nrh ; i++) {
134 retval[i] = retval[i-1] + width;
135 retval[i][ncl] = retval[i-1][ncl] + width*depth; /* first cell in column */
136 for(long j = ncl+1 ; j <= nch ; j++) {
137 retval[i][j] = retval[i][j-1] + depth;
138 }
139 assert(retval[i][nch]-retval[i][ncl] == (nch-ncl)*depth);
140 }
141 assert(retval[nrh]-retval[nrl] == (nrh-nrl)*width);
142 assert(&retval[nrh][nch][ndh]-&retval[nrl][ncl][ndl] == (nrh-nrl+1)*(nch-ncl+1)*(ndh-ndl+1)-1);
143
144 return retval;
145}
146
147/*--------------------------------------------------------------------------*/
148void
149free_ivector (int *v, long nl, long nh)
150/* free an int vector allocated with ivector() */
151{
152 free(v+nl);
153}
154
155/*--------------------------------------------------------------------------*/
156void
157free_dvector (double *v, long nl, long nh)
158/* free an double vector allocated with dvector() */
159{
160 free(v+nl);
161}
162
163/*--------------------------------------------------------------------------*/
164void
165free_imatrix (int **m, long nrl, long nrh, long ncl, long nch)
166/* free an int matrix allocated by imatrix() */
167{
168 free(m[nrl]+ncl);
169 free(m+nrl);
170}
171
172/*--------------------------------------------------------------------------*/
173void
174free_dmatrix (double **m, long nrl, long nrh, long ncl, long nch)
175/* free a double matrix allocated by dmatrix() */
176{
177 free(m[nrl]+ncl);
178 free(m+nrl);
179}
180
181/*--------------------------------------------------------------------------*/
182void
183free_d3tensor (double ***t, long nrl, long nrh, long ncl, long nch,
184 long ndl, long ndh)
185/* free a double d3tensor allocated by d3tensor() */
186{
187 free(t[nrl][ncl]+ndl);
188 free(t[nrl]+ncl);
189 free(t+nrl);
190}
191
192/*--------------------------------------------------------------------------*/
193int
194minimum2 (int i, int j)
195{
196 int result = i;
197 if (j < result)
198 result = j;
199 return result;
200}
201
202/*-------------------------------------------------------------------------*/
203int
204minimum3 (int i, int j, int k)
205{
206 int result = i;
207 if (j < result)
208 result = j;
209 if (k < result)
210 result = k;
211 return result;
212}
213
214/*--------------------------------------------------------------------------*/
215int
216maximum2 (int i, int j)
217{
218 int result = i;
219 if (j > result)
220 result = j;
221 return result;
222}
223
224/*--------------------------------------------------------------------------*/
225int
226maximum3 (int i, int j, int k)
227{
228 int result = i;
229 if (j > result)
230 result = j;
231 if (k > result)
232 result = k;
233 return result;
234}
235
236/*--------------------------------------------------------------------------*/
237int
238pow_int (int mantisse, int exponent)
239{
240 int i, result = 1;
241
242 for (i = 1; i <= exponent; i++)
243 result *= mantisse;
244
245 return result;
246}
247
248/*--------------------------------------------------------------------------*/
249void
250chebft_Zeros (double u[], int n, int inv)
251 /* eq. 5.8.7 and 5.8.8 at x = (5.8.4) of 2nd edition C++ NR */
252{
253 int k, j, isignum;
254 double fac, sum, Pion, *c;
255
256 c = dvector (0, n);
257 Pion = Pi / n;
258 if (inv == 0)
259 {
260 fac = 2.0 / n;
261 isignum = 1;
262 for (j = 0; j < n; j++)
263 {
264 sum = 0.0;
265 for (k = 0; k < n; k++)
266 sum += u[k] * cos (Pion * j * (k + 0.5));
267 c[j] = fac * sum * isignum;
268 isignum = -isignum;
269 }
270 }
271 else
272 {
273 for (j = 0; j < n; j++)
274 {
275 sum = -0.5 * u[0];
276 isignum = 1;
277 for (k = 0; k < n; k++)
278 {
279 sum += u[k] * cos (Pion * (j + 0.5) * k) * isignum;
280 isignum = -isignum;
281 }
282 c[j] = sum;
283 }
284 }
285 for (j = 0; j < n; j++)
286#if 0
287 if (fabs(c[j]) < 5.e-16)
288 u[j] = 0.0;
289 else
290#endif
291 u[j] = c[j];
292 free_dvector (c, 0, n);
293}
294
295/* --------------------------------------------------------------------------*/
296
297void
298chebft_Extremes (double u[], int n, int inv)
299 /* eq. 5.8.7 and 5.8.8 at x = (5.8.5) of 2nd edition C++ NR */
300{
301 int k, j, isignum, N = n - 1;
302 double fac, sum, PioN, *c;
303
304 c = dvector (0, N);
305 PioN = Pi / N;
306 if (inv == 0)
307 {
308 fac = 2.0 / N;
309 isignum = 1;
310 for (j = 0; j < n; j++)
311 {
312 sum = 0.5 * (u[0] + u[N] * isignum);
313 for (k = 1; k < N; k++)
314 sum += u[k] * cos (PioN * j * k);
315 c[j] = fac * sum * isignum;
316 isignum = -isignum;
317 }
318 c[N] = 0.5 * c[N];
319 }
320 else
321 {
322 for (j = 0; j < n; j++)
323 {
324 sum = -0.5 * u[0];
325 isignum = 1;
326 for (k = 0; k < n; k++)
327 {
328 sum += u[k] * cos (PioN * j * k) * isignum;
329 isignum = -isignum;
330 }
331 c[j] = sum;
332 }
333 }
334 for (j = 0; j < n; j++)
335 u[j] = c[j];
336 free_dvector (c, 0, N);
337}
338
339/* --------------------------------------------------------------------------*/
340
341void
342chder (double *c, double *cder, int n)
343{
344 int j;
345
346 cder[n] = 0.0;
347 cder[n - 1] = 0.0;
348 for (j = n - 2; j >= 0; j--)
349 cder[j] = cder[j + 2] + 2 * (j + 1) * c[j + 1];
350}
351
352/* --------------------------------------------------------------------------*/
353double
354chebev (double a, double b, double c[], int m, double x)
355 /* eq. 5.8.11 of C++ NR (2nd ed) */
356{
357 int j;
358 double djp2, djp1, dj; /* d_{j+2}, d_{j+1} and d_j */
359 double y;
360
361 /* rescale input to lie within [-1,1] */
362 y = 2*(x - 0.5*(b+a))/(b-a);
363
364 dj = djp1 = 0;
365 for(j = m-1 ; j >= 1; j--)
366 {
367 /* advance the coefficients */
368 djp2 = djp1;
369 djp1 = dj;
370 dj = 2*y*djp1 - djp2 + c[j];
371 }
372
373 return y*dj - djp1 + 0.5*c[0];
374}
375
376/* --------------------------------------------------------------------------*/
377void
378fourft (double *u, int N, int inv)
379 /* a (slow) Fourier transform, seems to be just eq. 12.1.6 and 12.1.9 of C++ NR (2nd ed) */
380{
381 int l, k, iy, M;
382 double x, x1, fac, Pi_fac, *a, *b;
383
384 M = N / 2;
385 a = dvector (0, M);
386 b = dvector (1, M); /* Actually: b=vector(1,M-1) but this is problematic if M=1*/
387 fac = 1. / M;
388 Pi_fac = Pi * fac;
389 if (inv == 0)
390 {
391 for (l = 0; l <= M; l++)
392 {
393 a[l] = 0;
394 if (l > 0 && l < M)
395 b[l] = 0;
396 x1 = Pi_fac * l;
397 for (k = 0; k < N; k++)
398 {
399 x = x1 * k;
400 a[l] += fac * u[k] * cos (x);
401 if (l > 0 && l < M)
402 b[l] += fac * u[k] * sin (x);
403 }
404 }
405 u[0] = a[0];
406 u[M] = a[M];
407 for (l = 1; l < M; l++)
408 {
409 u[l] = a[l];
410 u[l + M] = b[l];
411 }
412 }
413 else
414 {
415 a[0] = u[0];
416 a[M] = u[M];
417 for (l = 1; l < M; l++)
418 {
419 a[l] = u[l];
420 b[l] = u[M + l];
421 }
422 iy = 1;
423 for (k = 0; k < N; k++)
424 {
425 u[k] = 0.5 * (a[0] + a[M] * iy);
426 x1 = Pi_fac * k;
427 for (l = 1; l < M; l++)
428 {
429 x = x1 * l;
430 u[k] += a[l] * cos (x) + b[l] * sin (x);
431 }
432 iy = -iy;
433 }
434 }
435 free_dvector (a, 0, M);
436 free_dvector (b, 1, M);
437}
438
439/* -----------------------------------------*/
440void
441fourder (double u[], double du[], int N)
442{
443 int l, M, lpM;
444
445 M = N / 2;
446 du[0] = 0.;
447 du[M] = 0.;
448 for (l = 1; l < M; l++)
449 {
450 lpM = l + M;
451 du[l] = u[lpM] * l;
452 du[lpM] = -u[l] * l;
453 }
454}
455
456/* -----------------------------------------*/
457void
458fourder2 (double u[], double d2u[], int N)
459{
460 int l, l2, M, lpM;
461
462 d2u[0] = 0.;
463 M = N / 2;
464 for (l = 1; l <= M; l++)
465 {
466 l2 = l * l;
467 lpM = l + M;
468 d2u[l] = -u[l] * l2;
469 if (l < M)
470 d2u[lpM] = -u[lpM] * l2;
471 }
472}
473
474/* ----------------------------------------- */
475double
476fourev (double *u, int N, double x)
477{
478 int l, M = N / 2;
479 double xl, result;
480
481 result = 0.5 * (u[0] + u[M] * cos (x * M));
482 for (l = 1; l < M; l++)
483 {
484 xl = x * l;
485 result += u[l] * cos (xl) + u[M + l] * sin (xl);
486 }
487 return result;
488}
489
490/* ------------------------------------------------------------------------*/
491double
492norm1 (double *v, int n)
493{
494 int i;
495 double result = -1;
496
497 for (i = 0; i < n; i++)
498 if (fabs (v[i]) > result)
499 result = fabs (v[i]);
500
501 return result;
502}
503
504/* -------------------------------------------------------------------------*/
505double
506norm2 (double *v, int n)
507{
508 int i;
509 double result = 0;
510
511 for (i = 0; i < n; i++)
512 result += v[i] * v[i];
513
514 return sqrt (result);
515}
516
517/* -------------------------------------------------------------------------*/
518double
519scalarproduct (double *v, double *w, int n)
520{
521 int i;
522 double result = 0;
523
524 for (i = 0; i < n; i++)
525 result += v[i] * w[i];
526
527 return result;
528}
529
530/* -------------------------------------------------------------------------*/
531
532} // namespac Utilities
533} // namespace TP
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
#define Pi
Another nice thing: The TwoPunctures Error, borrowed from Pizza.
Definition TP_Logging.h:66
int minimum3(int i, int j, int k)
void fourft(double *u, int N, int inv)
double norm2(double *v, int n)
double chebev(double a, double b, double c[], int m, double x)
int maximum3(int i, int j, int k)
int maximum2(int i, int j)
double * dvector(long nl, long nh)
void chebft_Extremes(double u[], int n, int inv)
void free_dvector(double *v, long nl, long nh)
double norm1(double *v, int n)
double scalarproduct(double *v, double *w, int n)
void free_ivector(int *v, long nl, long nh)
void chebft_Zeros(double u[], int n, int inv)
void fourder2(double u[], double d2u[], int N)
int minimum2(int i, int j)
void chder(double *c, double *cder, int n)
double *** d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
int * ivector(long nl, long nh)
void fourder(double u[], double du[], int N)
double fourev(double *u, int N, double x)
double ** dmatrix(long nrl, long nrh, long ncl, long nch)
void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
int pow_int(int mantisse, int exponent)
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
int ** imatrix(long nrl, long nrh, long ncl, long nch)
This file contains aliases for making access to the long state vector Q as used eg.