aboutsummaryrefslogtreecommitdiff
path: root/src/interpolate.cc
blob: 5c9946c742809f5b25466e820eea974539de8d5d (plain)
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
#include <stdio.h>

#include "interpolate.hh"

static
void report_interp_error(int ierr)
{
  if (ierr < 0)
  {
    CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
        "CCTK_InterpGridArrays returned error code %d",ierr);
  }
}
  
void Multipole_Interp(CCTK_ARGUMENTS,
                      CCTK_REAL xs[], CCTK_REAL ys[], CCTK_REAL zs[],
                      int real_idx, int imag_idx,
                      CCTK_REAL sphere_real[], CCTK_REAL sphere_imag[])
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  // Need parameters for the following:
  // ntheta (dtheta = pi/(ntheta) 
  // nphi (dphi = 2pi/(nphi)
  // r (radius of sphere)
  // NOTE: depending on the interval of integration, denominator above may
  // need to be modified to avoid double counting


  CCTK_INT num_input_arrays  = imag_idx == -1 ? 1 : 2;
  CCTK_INT num_output_arrays = imag_idx == -1 ? 1 : 2;
  const CCTK_INT num_dims = 3;
  int ierr = -1;

  const void* interp_coords[num_dims] 
    = { (const void *) xs,
        (const void *) ys,
        (const void *) zs };

  const CCTK_INT input_array_indices[2]
    = { real_idx,
        imag_idx };

  const CCTK_INT output_array_types[2]
    = { CCTK_VARIABLE_REAL,
        CCTK_VARIABLE_REAL };

  void * output_arrays[2]
    = { (void *) sphere_real,
        (void *) sphere_imag };

  const int operator_handle = CCTK_InterpHandle(interpolator_name);

  int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
  ierr = Util_TableSetFromString(param_table_handle, interpolator_pars);

  const int coord_system_handle = CCTK_CoordSystemHandle(coord_system);

  ierr = CCTK_InterpGridArrays(
      cctkGH,
      num_dims,
      operator_handle,
      param_table_handle,
      coord_system_handle,
      CCTK_MyProc(cctkGH) == 0 ? (ntheta+1)*(nphi+1) : 0, // Only the 0 processor needs the points
      CCTK_VARIABLE_REAL,
      interp_coords,
      num_input_arrays,
      input_array_indices,
      num_output_arrays,
      output_array_types,
      output_arrays);

  report_interp_error(ierr);

  if (imag_idx == -1)
  {
    for (int i = 0; i < (ntheta+1)*(nphi+1); i++)
    {
      sphere_imag[i] = 0;
    }
  }

  Util_TableDestroy(param_table_handle);
}

// // Debugging routine
// void Multipole_InterpVar(CCTK_ARGUMENTS,
//                          CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[], const char *var_name,
//                          CCTK_REAL sphere_var[])
// {
//   DECLARE_CCTK_ARGUMENTS;
//   DECLARE_CCTK_PARAMETERS;

//   // Need parameters for the following:
//   // ntheta (dtheta = pi/(ntheta) 
//   // nphi (dphi = 2pi/(nphi)
//   // r (radius of sphere)
//   // NOTE: depending on the interval of integration, denominator above may
//   // need to be modified to avoid double counting

// //  printf("Multipole: Interpolating %s\n", var_name);

//   const CCTK_INT num_input_arrays = 1;
//   const CCTK_INT num_output_arrays = 1;
//   const CCTK_INT num_dims = 3;
//   int ierr = -1;

//   const void* interp_coords[num_dims] 
//     = { (const void *) x,
//         (const void *) y,
//         (const void *) z };

//   const CCTK_INT input_array_indices[num_input_arrays]
//     = { CCTK_VarIndex(var_name) };

//   const CCTK_INT output_array_types[num_output_arrays]
//     = { CCTK_VARIABLE_REAL };

//   void * output_arrays[num_output_arrays]
//     = { (void *) sphere_var };

//   const int operator_handle = CCTK_InterpHandle(interpolator_name);

//   int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
//   ierr = Util_TableSetFromString(param_table_handle, interpolator_pars);

//   const int coord_system_handle = CCTK_CoordSystemHandle(coord_system);

//   ierr = CCTK_InterpGridArrays(
//       cctkGH,
//       num_dims,
//       operator_handle,
//       param_table_handle,
//       coord_system_handle,
//       (ntheta+1)*(nphi+1),
//       CCTK_VARIABLE_REAL,
//       interp_coords,
//       num_input_arrays,
//       input_array_indices,
//       num_output_arrays,
//       output_array_types,
//       output_arrays);
//   report_interp_error(ierr);
// }