HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
IBCreateFacetMapping.c
Go to the documentation of this file.
1 
6 #include <stdlib.h>
7 #include <string.h>
8 #include <arrayfunctions.h>
9 #include <mathfunctions.h>
10 #include <immersedboundaries.h>
11 #include <mpivars.h>
12 
14 static inline int isInside(
15  double x,
16  double a,
17  double b
18  )
19 {
20  return ((x >= a) && (x <= b));
21 }
22 
28 static int interpNodesCoeffs(
29  void *m,
30  double xc,
31  double yc,
32  double zc,
33  double *x,
34  double *y,
35  double *z,
36  int *dim,
37  int ghosts,
38  char *mode,
39  int *ii,
41  int *jj,
43  int *kk,
45  int *inodes,
46  double *icoeffs
47  )
48 {
49  MPIVariables *mpi = (MPIVariables*) m;
50 
51  int i, j, k, ic, jc, kc;
52  ic = jc = kc = -1;
53 
54  double xmin = 0.5 * (x[ghosts-1] + x[ghosts]),
55  xmax = 0.5 * (x[dim[0]+ghosts-1] + x[dim[0]+ghosts]),
56  ymin = 0.5 * (y[ghosts-1] + y[ghosts]),
57  ymax = 0.5 * (y[dim[1]+ghosts-1] + y[dim[1]+ghosts]),
58  zmin = 0.5 * (z[ghosts-1] + z[ghosts]),
59  zmax = 0.5 * (z[dim[2]+ghosts-1] + z[dim[2]+ghosts]);
60 
61  for (i = 0; i < dim[0]+2*ghosts-1; i++) {
62  if (isInside(xc,x[i],x[i+1])) {
63  ic = i;
64  break;
65  }
66  }
67  if (ic <= ghosts-1) ic = ghosts;
68  else if (ic >= dim[0]+ghosts-1) ic = dim[0]+ghosts-2;
69 
70  for (j = 0; j < dim[1]+2*ghosts-1; j++) {
71  if (isInside(yc,y[j],y[j+1])) {
72  jc = j;
73  break;
74  }
75  }
76  if (jc <= ghosts-1) jc = ghosts;
77  else if (jc >= dim[1]+ghosts-1) jc = dim[1]+ghosts-2;
78 
79  for (k = 0; k < dim[2]+2*ghosts-1; k++) {
80  if (isInside(zc,z[k],z[k+1])) {
81  kc = k;
82  break;
83  }
84  }
85  if (kc <= ghosts-1) kc = ghosts;
86  else if (kc >= dim[2]+ghosts-1) kc = dim[2]+ghosts-2;
87 
88  if (!strcmp(mode,_IB_XY_)) { kc = ghosts; zc = 0.5*(zmin+zmax); }
89  else if (!strcmp(mode,_IB_XZ_)) { jc = ghosts; yc = 0.5*(ymin+ymax); }
90  else if (!strcmp(mode,_IB_YZ_)) { ic = ghosts; xc = 0.5*(xmin+xmax); }
91 
92  if (ic == -1) {
93  fprintf(stderr,"Error in interpNodesCoeffs() (in ImmersedBoundaries/IBCreateFacetMapping.c) on rank %d: ic = -1.\n", mpi->rank);
94  return(1);
95  }
96  if (jc == -1) {
97  fprintf(stderr,"Error in interpNodesCoeffs() (in ImmersedBoundaries/IBCreateFacetMapping.c) on rank %d: jc = -1.\n", mpi->rank);
98  return(1);
99  }
100  if (kc == -1) {
101  fprintf(stderr,"Error in interpNodesCoeffs() (in ImmersedBoundaries/IBCreateFacetMapping.c) on rank %d: kc = -1.\n", mpi->rank);
102  return(1);
103  }
104  ic++;
105  jc++;
106  kc++;
107 
108  if (ii) *ii = ic;
109  if (jj) *jj = jc;
110  if (kk) *kk = kc;
111 
112  int pc[_IB_NNODES_], index[_IB_NDIMS_];
113  index[0]=ic-1-ghosts; index[1]=jc-1-ghosts; index[2]=kc-1-ghosts; _ArrayIndex1D_(_IB_NDIMS_,dim,index,ghosts,pc[0]);
114  index[0]=ic-ghosts ; index[1]=jc-1-ghosts; index[2]=kc-1-ghosts; _ArrayIndex1D_(_IB_NDIMS_,dim,index,ghosts,pc[1]);
115  index[0]=ic-1-ghosts; index[1]=jc-ghosts ; index[2]=kc-1-ghosts; _ArrayIndex1D_(_IB_NDIMS_,dim,index,ghosts,pc[2]);
116  index[0]=ic-ghosts ; index[1]=jc-ghosts ; index[2]=kc-1-ghosts; _ArrayIndex1D_(_IB_NDIMS_,dim,index,ghosts,pc[3]);
117  index[0]=ic-1-ghosts; index[1]=jc-1-ghosts; index[2]=kc-ghosts ; _ArrayIndex1D_(_IB_NDIMS_,dim,index,ghosts,pc[4]);
118  index[0]=ic-ghosts ; index[1]=jc-1-ghosts; index[2]=kc-ghosts ; _ArrayIndex1D_(_IB_NDIMS_,dim,index,ghosts,pc[5]);
119  index[0]=ic-1-ghosts; index[1]=jc-ghosts ; index[2]=kc-ghosts ; _ArrayIndex1D_(_IB_NDIMS_,dim,index,ghosts,pc[6]);
120  index[0]=ic-ghosts ; index[1]=jc-ghosts ; index[2]=kc-ghosts ; _ArrayIndex1D_(_IB_NDIMS_,dim,index,ghosts,pc[7]);
121  _ArrayCopy1D_(pc,inodes,_IB_NNODES_);
122 
123  double coeffs[_IB_NNODES_];
124  TrilinearInterpCoeffs(x[ic-1],x[ic],y[jc-1],y[jc],z[kc-1],z[kc],xc,yc,zc,&coeffs[0]);
125  _ArrayCopy1D_(coeffs,icoeffs,_IB_NNODES_);
126 
127  return(0);
128 }
129 
143 int IBCreateFacetMapping(
144  void *ib,
145  void *m,
146  double *X,
147  int *dim,
148  int ghosts
149  )
150 {
151  ImmersedBoundary *IB = (ImmersedBoundary*) ib;
152  MPIVariables *mpi = (MPIVariables*) m;
153  Body3D *body = IB->body;
154  int nfacets = body->nfacets, n, count, ierr;
155  Facet3D *facets = body->surface;
156 
157  double *x = X,
158  *y = (x + dim[0] + 2*ghosts),
159  *z = (y + dim[1] + 2*ghosts);
160 
161  double xmin = 0.5 * (x[ghosts-1] + x[ghosts]),
162  xmax = 0.5 * (x[dim[0]+ghosts-1] + x[dim[0]+ghosts]),
163  ymin = 0.5 * (y[ghosts-1] + y[ghosts]),
164  ymax = 0.5 * (y[dim[1]+ghosts-1] + y[dim[1]+ghosts]),
165  zmin = 0.5 * (z[ghosts-1] + z[ghosts]),
166  zmax = 0.5 * (z[dim[2]+ghosts-1] + z[dim[2]+ghosts]);
167 
168  count = 0;
169  for (n = 0; n < nfacets; n++) {
170 
171  /* find facet centroid */
172  double xc, yc, zc;
173  xc = (facets[n].x1 + facets[n].x2 + facets[n].x3) / 3.0;
174  yc = (facets[n].y1 + facets[n].y2 + facets[n].y3) / 3.0;
175  zc = (facets[n].z1 + facets[n].z2 + facets[n].z3) / 3.0;
176 
177  if (!strcmp(IB->mode,_IB_3D_)) {
178  if (isInside(xc,xmin,xmax) && isInside(yc,ymin,ymax) && isInside(zc,zmin,zmax)) count++;
179  } else if (!strcmp(IB->mode,_IB_XY_)) {
180  if (isInside(xc,xmin,xmax) && isInside(yc,ymin,ymax)) count++;
181  } else if (!strcmp(IB->mode,_IB_XZ_)) {
182  if (isInside(xc,xmin,xmax) && isInside(zc,zmin,zmax)) count++;
183  } else if (!strcmp(IB->mode,_IB_YZ_)) {
184  if (isInside(yc,ymin,ymax) && isInside(zc,zmin,zmax)) count++;
185  }
186  }
187 
188  int nfacets_local = count;
189  if (nfacets_local > 0) {
190 
191  FacetMap *fmap = (FacetMap*) calloc (nfacets_local, sizeof(FacetMap));
192  count = 0;
193 
194  for (n = 0; n < nfacets; n++) {
195 
196  double xc, yc, zc;
197  xc = (facets[n].x1 + facets[n].x2 + facets[n].x3) / 3.0;
198  yc = (facets[n].y1 + facets[n].y2 + facets[n].y3) / 3.0;
199  zc = (facets[n].z1 + facets[n].z2 + facets[n].z3) / 3.0;
200 
201  int flag = 0;
202  if (!strcmp(IB->mode,_IB_3D_)) {
203  if (isInside(xc,xmin,xmax) && isInside(yc,ymin,ymax) && isInside(zc,zmin,zmax)) flag = 1;
204  } else if (!strcmp(IB->mode,_IB_XY_)) {
205  if (isInside(xc,xmin,xmax) && isInside(yc,ymin,ymax)) flag = 1;
206  } else if (!strcmp(IB->mode,_IB_XZ_)) {
207  if (isInside(xc,xmin,xmax) && isInside(zc,zmin,zmax)) flag = 1;
208  } else if (!strcmp(IB->mode,_IB_YZ_)) {
209  if (isInside(yc,ymin,ymax) && isInside(zc,zmin,zmax)) flag = 1;
210  }
211 
212  if (flag == 1) {
213  fmap[count].facet = facets + n;
214  fmap[count].index = n;
215 
216  fmap[count].xc = xc;
217  fmap[count].yc = yc;
218  fmap[count].zc = zc;
219 
220  int ic,jc, kc;
221 
222  ierr = interpNodesCoeffs( mpi,
223  xc, yc, zc,
224  x, y, z,
225  dim,
226  ghosts,
227  IB->mode,
228  &ic, &jc, &kc,
229  fmap[count].interp_nodes,
230  fmap[count].interp_coeffs );
231  if (ierr) {
232  fprintf(stderr, "Error in IBCreateFacetMapping(): \n");
233  fprintf(stderr, " interpNodesCoeffs() returned with error code %d on rank %d.\n",
234  ierr, mpi->rank );
235  return(ierr);
236  }
237 
238  double dx = x[ic] - x[ic-1];
239  double dy = y[jc] - y[jc-1];
240  double dz = z[kc] - z[kc-1];
241  double ds;
242  if (!strcmp(IB->mode,_IB_XY_)) ds = min(dx,dy);
243  else if (!strcmp(IB->mode,_IB_XZ_)) ds = min(dx,dz);
244  else if (!strcmp(IB->mode,_IB_YZ_)) ds = min(dy,dz);
245  else ds = min3(dx,dy,dz);
246 
247  double nx = fmap[count].facet->nx;
248  double ny = fmap[count].facet->ny;
249  double nz = fmap[count].facet->nz;
250 
251  if (nx == 0.0) nx += IB->delta*ds;
252  if (ny == 0.0) ny += IB->delta*ds;
253  if (nz == 0.0) nz += IB->delta*ds;
254 
255  double xns = xc + sign(nx)*ds;
256  double yns = yc + sign(ny)*ds;
257  double zns = zc + sign(nz)*ds;
258 
259  fmap[count].dx = xns - xc;
260  fmap[count].dy = yns - yc;
261  fmap[count].dz = zns - zc;
262 
263  ierr = interpNodesCoeffs( mpi,
264  xns, yns, zns,
265  x, y, z,
266  dim,
267  ghosts,
268  IB->mode,
269  NULL,NULL,NULL,
270  fmap[count].interp_nodes_ns,
271  fmap[count].interp_coeffs_ns );
272  if (ierr) {
273  fprintf(stderr, "Error in IBCreateFacetMapping(): \n");
274  fprintf(stderr, " interpNodesCoeffs() returned with error code %d on rank %d.\n",
275  ierr, mpi->rank );
276  return(ierr);
277  }
278 
279  count++;
280  }
281  }
282 
283  IB->nfacets_local = nfacets_local;
284  IB->fmap = fmap;
285 
286  } else {
287 
288  IB->nfacets_local = 0;
289  IB->fmap = NULL;
290 
291  }
292 
293  return(0);
294 }
FacetMap::interp_coeffs
double interp_coeffs[_IB_NNODES_]
Definition: immersedboundaries.h:71
isInside
static int isInside(double x, double a, double b)
Definition: IBCreateFacetMapping.c:14
Facet3D
Structure defining a facet.
Definition: immersedboundaries.h:39
FacetMap
Structure defining a facet map.
Definition: immersedboundaries.h:67
Body3D
Structure defining a body.
Definition: immersedboundaries.h:100
_IB_NNODES_
#define _IB_NNODES_
Definition: immersedboundaries.h:13
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
mpivars.h
MPI related function definitions.
mathfunctions.h
Contains function definitions for common mathematical functions.
min
#define min(a, b)
Definition: math_ops.h:14
FacetMap::interp_coeffs_ns
double interp_coeffs_ns[_IB_NNODES_]
Definition: immersedboundaries.h:73
Body3D::surface
Facet3D * surface
Definition: immersedboundaries.h:102
_IB_3D_
#define _IB_3D_
Definition: immersedboundaries.h:22
FacetMap::interp_nodes
int interp_nodes[_IB_NNODES_]
Definition: immersedboundaries.h:70
sign
#define sign(a)
Definition: math_ops.h:54
TrilinearInterpCoeffs
void TrilinearInterpCoeffs(double, double, double, double, double, double, double, double, double, double *)
Definition: TrilinearInterpolation.c:22
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
FacetMap::facet
Facet3D * facet
Definition: immersedboundaries.h:68
_IB_NDIMS_
#define _IB_NDIMS_
Definition: immersedboundaries.h:10
_IB_YZ_
#define _IB_YZ_
Definition: immersedboundaries.h:20
immersedboundaries.h
Structures and function definitions for immersed boundaries.
interpNodesCoeffs
static int interpNodesCoeffs(void *m, double xc, double yc, double zc, double *x, double *y, double *z, int *dim, int ghosts, char *mode, int *ii, int *jj, int *kk, int *inodes, double *icoeffs)
Definition: IBCreateFacetMapping.c:28
_IB_XZ_
#define _IB_XZ_
Definition: immersedboundaries.h:18
FacetMap::interp_nodes_ns
int interp_nodes_ns[_IB_NNODES_]
Definition: immersedboundaries.h:72
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
ImmersedBoundary::mode
char mode[_MAX_STRING_SIZE_]
Definition: immersedboundaries.h:164
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
arrayfunctions.h
Contains macros and function definitions for common array operations.
min3
#define min3(a, b, c)
Definition: math_ops.h:23
IBCreateFacetMapping
int IBCreateFacetMapping(void *ib, void *m, double *X, int *dim, int ghosts)
Definition: IBCreateFacetMapping.c:143
_IB_XY_
#define _IB_XY_
Definition: immersedboundaries.h:16