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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
|
/*@@
@file Trace.c
@date Thu Apr 25 16:35:04 2002
@author Tom Goodale
@desc
Calculates the trace of things
@enddesc
@version $Header$
@@*/
#include "cctk.h"
#include "ADMAnalysis.h"
static const char *rcsid = "$Header$";
CCTK_FILEVERSION(CactusEinstein_ADMAnalysis_Trace_c)
/********************************************************************
********************* Local Data Types ***********************
********************************************************************/
#define SQR(a) ((a)*(a))
/********************************************************************
********************* Local Routine Prototypes *********************
********************************************************************/
/********************************************************************
***************** Scheduled Routine Prototypes *********************
********************************************************************/
/********************************************************************
********************* Other Routine Prototypes *********************
********************************************************************/
/********************************************************************
********************* Local Data *****************************
********************************************************************/
/********************************************************************
********************* External Routines **********************
********************************************************************/
/*@@
@routine ADMAnalysis_Trace
@date Thu Apr 25 16:49:57 2002
@author Tom Goodale
@desc
Calculates the trace of a symmetric 2 index 3-tensor.
Optionally returns the trace of the metric at each point.
@enddesc
@calls
@calledby
@history
@hdate Thu Apr 25 16:50:38 2002 @hauthor Tom Goodale
@hdesc Took old evlatrK routine and made general
@endhistory
@var lsh
@vdesc grid size
@vtype const int *
@vio in
@var g11
@vdesc the 11 component of the metric tensor
@vtype const CCTK_REAL *
@vio in
@var g12
@vdesc the 12 component of the metric tensor
@vtype const CCTK_REAL *
@vio in
@var g13
@vdesc the 13 component of the metric tensor
@vtype const CCTK_REAL *
@vio in
@var g22
@vdesc the 22 component of the metric tensor
@vtype const CCTK_REAL *
@vio in
@var g23
@vdesc the 23 component of the metric tensor
@vtype const CCTK_REAL *
@vio in
@var g33
@vdesc the 33 component of the metric tensor
@vtype const CCTK_REAL *
@vio in
@var tensor11
@vdesc the 11 component of the tensor
@vtype CCTK_REAL *
@vio out
@var tensor12
@vdesc the 12 component of the tensor
@vtype CCTK_REAL *
@vio out
@var tensor13
@vdesc the 13 component of the tensor
@vtype CCTK_REAL *
@vio out
@var tensor22
@vdesc the 22 component of the tensor
@vtype CCTK_REAL *
@vio out
@var tensor23
@vdesc the 23 component of the tensor
@vtype CCTK_REAL *
@vio out
@var tensor33
@vdesc the 33 component of the tensor
@vtype CCTK_REAL *
@vio out
@var detg
@vdesc the determinant of the metric tensor
@vtype CCTK_REAL *
@vio out
@vcomment
Will be ignored if NULL.
@endvar
@@*/
void ADMAnalysis_Trace(const CCTK_INT *lsh,
const CCTK_REAL *g11,
const CCTK_REAL *g12,
const CCTK_REAL *g13,
const CCTK_REAL *g22,
const CCTK_REAL *g23,
const CCTK_REAL *g33,
CCTK_REAL *tensor11,
CCTK_REAL *tensor12,
CCTK_REAL *tensor13,
CCTK_REAL *tensor22,
CCTK_REAL *tensor23,
CCTK_REAL *tensor33,
CCTK_REAL *trace,
CCTK_REAL *detg)
{
int i;
CCTK_REAL det, u11, u12, u22, u13, u23, u33, two;
CCTK_REAL g11p, g12p, g13p, g22p, g23p, g33p;
two = 2.0;
/* loop over all the gridpoints */
for(i = 0; i< lsh[0]*lsh[1]*lsh[2];i++)
{
/* get the metric */
g11p = g11[i];
g12p = g12[i];
g13p = g13[i];
g22p = g22[i];
g23p = g23[i];
g33p = g33[i];
/* compute determinant */
det=-(g13p*g13p*g22p)+2.*g12p*g13p*g23p-g11p*g23p*
g23p-g12p*g12p*g33p+g11p*g22p*g33p;
/* invert metric. This is the conformal upper metric */
u11=(-SQR(g23p) + g22p*g33p)/det;
u12=(g13p*g23p - g12p*g33p)/det;
u22=(-SQR(g13p) + g11p*g33p)/det;
u13=(-g13p*g22p + g12p*g23p)/det;
u23=(g12p*g13p - g11p*g23p)/det;
u33=(-SQR(g12p) + g11p*g22p)/det;
/* Calculate trK */
trace[i] = (u11*tensor11[i] + u22*tensor22[i] +
u33*tensor33[i]+ two*u12*tensor12[i] +
two*u13*tensor13[i] + two*u23*tensor23[i]);
if(detg)
{
detg[i]= det;
}
}
}
/********************************************************************
********************* Local Routines *************************
********************************************************************/
|
