HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
IBInterpCoeffs.c File Reference

Compute interpolation nodes and coefficients for immersed boundary points. More...

#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <mpivars.h>
#include <immersedboundaries.h>

Go to the source code of this file.

Functions

int IBInterpCoeffs (void *ib, void *m, double *X, int *dim_l, int ghosts, double *blank)
 

Detailed Description

Compute interpolation nodes and coefficients for immersed boundary points.

Author
Debojyoti Ghosh

Definition in file IBInterpCoeffs.c.

Function Documentation

int IBInterpCoeffs ( void *  ib,
void *  m,
double *  X,
int *  dim_l,
int  ghosts,
double *  blank 
)

Compute the interpolation nodes and coefficients for immersed boundary points: For each immersed boundary point, do the following:

  • From the immersed boundary point, extend a probe in the direction defined by the outward normal of the "nearest" facet (computed in IBNearestFacetNormal()), till the probe tip reaches a point in space such that all surrounding (_IB_NNODES_) grid points are "interior" points, i.e., outside the immersed body (they are "interior" to the computational domain).
  • Store the indices of the surrounding grid points, as well as the trilinear interpolation coefficients to interpolate a variable from the surrounding points to the probe tip.
Parameters
ibImmersed boundary object of type ImmersedBoundary
mMPI object of type MPIVariables
XArray of (local) spatial coordinates
dim_lInteger array of local grid size in each spatial dimension
ghostsNumber of ghost points
blankBlanking array: for grid points within the body, this value will be set to 0

Definition at line 22 of file IBInterpCoeffs.c.

31 {
32  ImmersedBoundary *IB = (ImmersedBoundary*) ib;
33  MPIVariables *mpi = (MPIVariables*) m;
34  IBNode *boundary = IB->boundary;
35 
36  double eps = IB->tolerance;
37  int maxiter = IB->itr_max,
38  n_boundary = IB->n_boundary_nodes;
39 
40  int imax = dim_l[0],
41  jmax = dim_l[1],
42  kmax = dim_l[2];
43 
44  int is = mpi->is[0],
45  js = mpi->is[1],
46  ks = mpi->is[2];
47 
48  static int index[_IB_NDIMS_];
49 
50  int dg;
51  for (dg = 0; dg < n_boundary; dg++) {
52  int i, j, k, p;
53  double xb, yb, zb;
54  double nx, ny, nz;
55  double xx, yy, zz;
56  double dx, dy, dz;
57  double ds, dist;
58  double xtip, ytip, ztip;
59 
60  index[0] = i = boundary[dg].i;
61  index[1] = j = boundary[dg].j;
62  index[2] = k = boundary[dg].k;
63  _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,p);
64 
65  xb = boundary[dg].x;
66  yb = boundary[dg].y;
67  zb = boundary[dg].z;
68 
69  nx = boundary[dg].face->nx;
70  ny = boundary[dg].face->ny;
71  nz = boundary[dg].face->nz;
72  xx = boundary[dg].face->x1;
73  yy = boundary[dg].face->y1;
74  zz = boundary[dg].face->z1;
75 
76  dist = nx*(xx-xb) + ny*(yy-yb) + nz*(zz-zb);
77 
78  double x1, x2, y1, y2, z1, z2;
79  _GetCoordinate_(0,(i+1),dim_l,ghosts,X,x1);
80  _GetCoordinate_(0,(i-1),dim_l,ghosts,X,x2);
81  _GetCoordinate_(1,(j+1),dim_l,ghosts,X,y1);
82  _GetCoordinate_(1,(j-1),dim_l,ghosts,X,y2);
83  _GetCoordinate_(2,(k+1),dim_l,ghosts,X,z1);
84  _GetCoordinate_(2,(k-1),dim_l,ghosts,X,z2);
85  dx = 0.5 * (x1 - x2);
86  dy = 0.5 * (y1 - y2);
87  dz = 0.5 * (z1 - z2);
88  ds = min3(dx, dy, dz);
89 
90  xtip = xb + dist*nx;
91  ytip = yb + dist*ny;
92  ztip = zb + dist*nz;
93 
94  int is_it_in = 0;
95  int iter = 0;
96  int itip, jtip, ktip;
97  while(!is_it_in && (iter < maxiter)) {
98  iter++;
99  itip = i;
100  jtip = j;
101  ktip = k;
102 
103  if (xtip > xb) {
104  double xx;
105  _GetCoordinate_(0,itip,dim_l,ghosts,X,xx);
106  while ((xx < xtip) && (itip < imax+ghosts-1)) {
107  itip++;
108  _GetCoordinate_(0,itip,dim_l,ghosts,X,xx);
109  }
110  } else {
111  double xx;
112  _GetCoordinate_(0,(itip-1),dim_l,ghosts,X,xx);
113  while ((xx > xtip) && (itip > -ghosts)) {
114  itip--;
115  _GetCoordinate_(0,(itip-1),dim_l,ghosts,X,xx);
116  }
117  }
118 
119  if (ytip > yb) {
120  double yy;
121  _GetCoordinate_(1,jtip,dim_l,ghosts,X,yy);
122  while ((yy < ytip) && (jtip < jmax+ghosts-1)) {
123  jtip++;
124  _GetCoordinate_(1,jtip,dim_l,ghosts,X,yy);
125  }
126  } else {
127  double yy;
128  _GetCoordinate_(1,(jtip-1),dim_l,ghosts,X,yy);
129  while ((yy > ytip) && (jtip > -ghosts)) {
130  jtip--;
131  _GetCoordinate_(1,(jtip-1),dim_l,ghosts,X,yy);
132  }
133  }
134 
135  if (ztip > zb) {
136  double zz;
137  _GetCoordinate_(2,ktip,dim_l,ghosts,X,zz);
138  while ((zz < ztip) && (ktip < kmax+ghosts-1)) {
139  ktip++;
140  _GetCoordinate_(2,ktip,dim_l,ghosts,X,zz);
141  }
142  } else {
143  double zz;
144  _GetCoordinate_(2,(ktip-1),dim_l,ghosts,X,zz);
145  while ((zz > ztip) && (ktip > -ghosts)) {
146  ktip--;
147  _GetCoordinate_(2,(ktip-1),dim_l,ghosts,X,zz);
148  }
149  }
150 
151  int ptip[_IB_NNODES_];
152  index[0] = itip ; index[1] = jtip ; index[2] = ktip ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[0]);
153  index[0] = itip-1; index[1] = jtip ; index[2] = ktip ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[1]);
154  index[0] = itip ; index[1] = jtip-1; index[2] = ktip ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[2]);
155  index[0] = itip ; index[1] = jtip ; index[2] = ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[3]);
156  index[0] = itip-1; index[1] = jtip-1; index[2] = ktip ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[4]);
157  index[0] = itip ; index[1] = jtip-1; index[2] = ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[5]);
158  index[0] = itip-1; index[1] = jtip ; index[2] = ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[6]);
159  index[0] = itip-1; index[1] = jtip-1; index[2] = ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[7]);
160 
161  int nflow = 0;
162  nflow += blank[ptip[0]];
163  nflow += blank[ptip[1]];
164  nflow += blank[ptip[2]];
165  nflow += blank[ptip[3]];
166  nflow += blank[ptip[4]];
167  nflow += blank[ptip[5]];
168  nflow += blank[ptip[6]];
169  nflow += blank[ptip[7]];
170  if (nflow == _IB_NNODES_) {
171  is_it_in = 1;
172  } else if (nflow < _IB_NNODES_) {
173  is_it_in = 0;
174  xtip += nx*absolute(ds);
175  ytip += ny*absolute(ds);
176  ztip += nz*absolute(ds);
177  } else {
178  fprintf(stderr,"Error in IBInterpCoeffs() (Bug in code) - counting interior points surrounding probe tip \n");
179  fprintf(stderr,"on rank %d.\n", mpi->rank);
180  fprintf(stderr,"Value of nflow is %d but can only be positive and <= %d.\n",nflow,_IB_NNODES_);
181  return(1);
182  }
183  }
184 
185  if (!is_it_in) {
186  fprintf(stderr,"Error in IBInterpCoeffs() on rank %d - interior point not found for immersed boundary point (%d,%d,%d)!\n",
187  mpi->rank, i, j, k);
188  return(1);
189  }
190 
191  double tlx[2],tly[2],tlz[2];
192  _GetCoordinate_(0,(itip-1),dim_l,ghosts,X,tlx[0]);
193  _GetCoordinate_(0,(itip ),dim_l,ghosts,X,tlx[1]);
194  _GetCoordinate_(1,(jtip-1),dim_l,ghosts,X,tly[0]);
195  _GetCoordinate_(1,(jtip ),dim_l,ghosts,X,tly[1]);
196  _GetCoordinate_(2,(ktip-1),dim_l,ghosts,X,tlz[0]);
197  _GetCoordinate_(2,(ktip ),dim_l,ghosts,X,tlz[1]);
198 
199  int ptip[_IB_NNODES_];
200  index[0]=itip-1; index[1]=jtip-1; index[2]=ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[0]);
201  index[0]=itip ; index[1]=jtip-1; index[2]=ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[1]);
202  index[0]=itip-1; index[1]=jtip ; index[2]=ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[2]);
203  index[0]=itip ; index[1]=jtip ; index[2]=ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[3]);
204  index[0]=itip-1; index[1]=jtip-1; index[2]=ktip ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[4]);
205  index[0]=itip ; index[1]=jtip-1; index[2]=ktip ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[5]);
206  index[0]=itip-1; index[1]=jtip ; index[2]=ktip ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[6]);
207  index[0]=itip ; index[1]=jtip ; index[2]=ktip ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[7]);
208  _ArrayCopy1D_(ptip,boundary[dg].interp_nodes,_IB_NNODES_);
209 
210  double coeffs[_IB_NNODES_];
211  TrilinearInterpCoeffs(tlx[0],tlx[1],tly[0],tly[1],tlz[0],tlz[1],xtip,ytip,ztip,&coeffs[0]);
212  _ArrayCopy1D_(coeffs,boundary[dg].interp_coeffs,_IB_NNODES_);
213 
214  double tipdist = absolute(nx*(xx-xtip) + ny*(yy-ytip) + nz*(zz-ztip));
215  boundary[dg].interp_node_distance = tipdist;
216  boundary[dg].surface_distance = absolute(dist);
217  if (tipdist < eps) {
218  fprintf(stderr,"Warning in IBInterpCoeffs() on rank %d - how can probe tip be on surface? Tipdist = %e\n",
219  mpi->rank,tipdist);
220  }
221 
222  }
223  return(0);
224 }
IBNode
Structure defining an immersed boundary node.
Definition: immersedboundaries.h:125
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
min3
#define min3(a, b, c)
Definition: math_ops.h:23
_IB_NDIMS_
#define _IB_NDIMS_
Definition: immersedboundaries.h:10
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
TrilinearInterpCoeffs
void TrilinearInterpCoeffs(double, double, double, double, double, double, double, double, double, double *)
Definition: TrilinearInterpolation.c:22
IBNode::interp_node_distance
double interp_node_distance
Definition: immersedboundaries.h:136
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
absolute
#define absolute(a)
Definition: math_ops.h:32
IBNode::surface_distance
double surface_distance
Definition: immersedboundaries.h:136
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
_IB_NNODES_
#define _IB_NNODES_
Definition: immersedboundaries.h:13
IBNode::face
Facet3D * face
Definition: immersedboundaries.h:133