HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
IBInterpCoeffs.c
Go to the documentation of this file.
1 
6 #include <basic.h>
7 #include <arrayfunctions.h>
8 #include <mathfunctions.h>
9 #include <mpivars.h>
10 #include <immersedboundaries.h>
11 
22 int IBInterpCoeffs(
23  void *ib,
24  void *m,
25  double *X,
26  int *dim_l,
27  int ghosts,
28  double *blank
30  )
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 }
absolute
#define absolute(a)
Definition: math_ops.h:32
_IB_NNODES_
#define _IB_NNODES_
Definition: immersedboundaries.h:13
IBNode
Structure defining an immersed boundary node.
Definition: immersedboundaries.h:125
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.
basic.h
Some basic definitions and macros.
IBNode::interp_node_distance
double interp_node_distance
Definition: immersedboundaries.h:136
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
_IB_NDIMS_
#define _IB_NDIMS_
Definition: immersedboundaries.h:10
IBNode::face
Facet3D * face
Definition: immersedboundaries.h:133
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
immersedboundaries.h
Structures and function definitions for immersed boundaries.
IBInterpCoeffs
int IBInterpCoeffs(void *ib, void *m, double *X, int *dim_l, int ghosts, double *blank)
Definition: IBInterpCoeffs.c:22
IBNode::surface_distance
double surface_distance
Definition: immersedboundaries.h:136
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
_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