diff options
Diffstat (limited to 'src/Trace.c')
| -rw-r--r-- | src/Trace.c | 180 |
1 files changed, 180 insertions, 0 deletions
diff --git a/src/Trace.c b/src/Trace.c new file mode 100644 index 0000000..241900b --- /dev/null +++ b/src/Trace.c @@ -0,0 +1,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 ************************* + ********************************************************************/ + |