HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
FillGhostCells.c
Go to the documentation of this file.
1 
6 #include <arrayfunctions.h>
7 
17 void fillGhostCells( const int* const a_dim,
18  const int a_ngpt,
19  double* const a_u,
20  const int a_nvars,
21  const int a_ndims,
22  const int* const a_periodic
23  )
24 {
25  for (int d = 0; d < a_ndims; d++) {
26 
27  int bounds[a_ndims];
28  _ArrayCopy1D_(a_dim, bounds, a_ndims);
29  bounds[d] = a_ngpt;
30 
31  int index[a_ndims];
32  _ArraySetValue_(index, a_ndims, 0);
33 
34  if (a_periodic[d]) {
35 
36  /* periodic case */
37 
38  int done = 0;
39  while (!done) {
40 
41  {
42  /* low end - face = 1 */
43 
44  int p_gpt = 0,
45  p_int = 0;
46 
47  int index_gpt[a_ndims];
48  _ArrayCopy1D_(index, index_gpt, a_ndims);
49  index_gpt[d] -= a_ngpt;
50  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
51 
52  int index_int[a_ndims];
53  _ArrayCopy1D_(index, index_int, a_ndims);
54  index_int[d] += (a_dim[d]-a_ngpt);
55  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int);
56 
57  _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
58  }
59 
60  {
61  /* high end - face = -1 */
62 
63  int p_gpt = 0,
64  p_int = 0;
65 
66  int index_gpt[a_ndims];
67  _ArrayCopy1D_(index, index_gpt, a_ndims);
68  index_gpt[d] += a_dim[d];
69  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
70 
71  int index_int[a_ndims];
72  _ArrayCopy1D_(index, index_int, a_ndims);
73  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int);
74 
75  _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
76  }
77 
78  _ArrayIncrementIndex_(a_ndims, bounds, index, done);
79 
80  }
81 
82  } else {
83 
84  /* not periodic - extrapolate */
85 
86  int done = 0;
87  while (!done) {
88 
89  {
90  /* low end - face = 1 */
91 
92  int p_gpt = 0,
93  p_int_0 = 0,
94  p_int_1 = 0,
95  p_int_2 = 0,
96  p_int_3 = 0;
97 
98  int index_gpt[a_ndims];
99  _ArrayCopy1D_(index, index_gpt, a_ndims);
100  index_gpt[d] -= a_ngpt;
101  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
102 
103  int index_int[a_ndims];
104  _ArrayCopy1D_(index, index_int, a_ndims);
105 
106  index_int[d] = 0;
107  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_0);
108  index_int[d]++;
109  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_1);
110  index_int[d]++;
111  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_2);
112  index_int[d]++;
113  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_3);
114 
115  double alpha = - (double) (a_ngpt - index[d]);
116  double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
117  double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
118  double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
119  double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;
120 
121  for (int v = 0; v < a_nvars; v++) {
122 
123  a_u[p_gpt*a_nvars+v] = c0 * a_u[p_int_0*a_nvars+v]
124  + c1 * a_u[p_int_1*a_nvars+v]
125  + c2 * a_u[p_int_2*a_nvars+v]
126  + c3 * a_u[p_int_3*a_nvars+v];
127 
128  }
129 
130  }
131 
132  {
133  /* high end - face = -1 */
134 
135  int p_gpt = 0,
136  p_int_0 = 0,
137  p_int_1 = 0,
138  p_int_2 = 0,
139  p_int_3 = 0;
140 
141  int index_gpt[a_ndims];
142  _ArrayCopy1D_(index, index_gpt, a_ndims);
143  index_gpt[d] += a_dim[d];
144  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
145 
146  int index_int[a_ndims];
147  _ArrayCopy1D_(index, index_int, a_ndims);
148 
149  index_int[d] = a_dim[d]-1;
150  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_0);
151  index_int[d]--;
152  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_1);
153  index_int[d]--;
154  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_2);
155  index_int[d]--;
156  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_3);
157 
158  double alpha = - (double) (index[d]+1);
159  double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
160  double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
161  double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
162  double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;
163 
164  for (int v = 0; v < a_nvars; v++) {
165 
166  a_u[p_gpt*a_nvars+v] = c0 * a_u[p_int_0*a_nvars+v]
167  + c1 * a_u[p_int_1*a_nvars+v]
168  + c2 * a_u[p_int_2*a_nvars+v]
169  + c3 * a_u[p_int_3*a_nvars+v];
170 
171  }
172 
173  }
174 
175  _ArrayIncrementIndex_(a_ndims, bounds, index, done);
176 
177  }
178 
179  }
180 
181  }
182 
183  return;
184 }
fillGhostCells
void fillGhostCells(const int *const a_dim, const int a_ngpt, double *const a_u, const int a_nvars, const int a_ndims, const int *const a_periodic)
Definition: FillGhostCells.c:17
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
arrayfunctions.h
Contains macros and function definitions for common array operations.