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
IBIdentifyBody.c File Reference

Identify grid points inside immersed body. More...

#include <stdio.h>
#include <stdlib.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <immersedboundaries.h>
#include <mpivars.h>

Go to the source code of this file.

Functions

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

Detailed Description

Identify grid points inside immersed body.

Author
Debojyoti Ghosh

Definition in file IBIdentifyBody.c.

Function Documentation

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

Identify the grid points of the given grid that are inside a given body whose surface is defined as an unstructured triangulation. This function uses the ray-tracing method and is a copy of a FORTRAN function originally written by Dr. Jay Sitaraman.

Parameters
ibImmersed boundary object of type ImmersedBoundary
dim_gglobal dimensions
dim_llocal dimensions
ghostsnumber of ghost points
mMPI object of type MPIVariables
XArray of global spatial coordinates
blankBlanking array: for grid points within the body, this value will be set to 0

Definition at line 20 of file IBIdentifyBody.c.

30 {
31  ImmersedBoundary *IB = (ImmersedBoundary*) ib;
32  MPIVariables *mpi = (MPIVariables*) m;
33  Body3D *body = IB->body;
34  double *x = X,
35  *y = X+dim_g[0],
36  *z = X+dim_g[0]+dim_g[1],
37  eps = IB->tolerance;
38  int itr_max = IB->itr_max,
39  i, j, k, n, v;
40 
41  double xmax = body->xmax,
42  xmin = body->xmin,
43  ymax = body->ymax,
44  ymin = body->ymin,
45  zmax = body->zmax,
46  zmin = body->zmin;
47 
48  double fac = 1.5;
49  double Lx, Ly, Lz;
50  Lx = (xmax - xmin);
51  Ly = (ymax - ymin);
52  Lz = (zmax - zmin);
53  double xc, yc, zc;
54  xc = 0.5 * (xmin + xmax);
55  yc = 0.5 * (ymin + ymax);
56  zc = 0.5 * (zmin + zmax);
57  xmax = xc + fac * Lx/2;
58  xmin = xc - fac * Lx/2;
59  ymax = yc + fac * Ly/2;
60  ymin = yc - fac * Ly/2;
61  zmax = zc + fac * Lz/2;
62  zmin = zc - fac * Lz/2;
63 
64  int imin, imax, jmin, jmax, kmin, kmax;
65  imin = dim_g[0]-1; imax = 0;
66  jmin = dim_g[1]-1; jmax = 0;
67  kmin = dim_g[2]-1; kmax = 0;
68  for (i = 0; i < dim_g[0]; i++) {
69  for (j = 0; j < dim_g[1]; j++) {
70  for (k = 0; k < dim_g[2]; k++) {
71  if ( ((x[i]-xmin)*(x[i]-xmax) < 0)
72  && ((y[j]-ymin)*(y[j]-ymax) < 0)
73  && ((z[k]-zmin)*(z[k]-zmax) < 0)) {
74  imin = min(i, imin);
75  imax = max(i, imax);
76  jmin = min(j, jmin);
77  jmax = max(j, jmax);
78  kmin = min(k, kmin);
79  kmax = max(k, kmax);
80  }
81  }
82  }
83  }
84 
85  double *cof[5], *xd, *dist;
86  xd = (double*) calloc (itr_max,sizeof(double));
87  dist = (double*) calloc (itr_max,sizeof(double));
88  for (v = 0; v < 5; v++) cof[v] = (double*) calloc(body->nfacets, sizeof(double));
89 
90  for (n = 0; n < body->nfacets; n++) {
91  cof[0][n] = body->surface[n].y1 - body->surface[n].y3;
92  cof[1][n] = body->surface[n].y2 - body->surface[n].y3;
93  cof[2][n] = body->surface[n].z1 - body->surface[n].z3;
94  cof[3][n] = body->surface[n].z2 - body->surface[n].z3;
95  double den = cof[0][n]*cof[3][n] - cof[1][n]*cof[2][n];
96  if (absolute(den) > eps) cof[4][n] = 1.0 / den;
97  else cof[4][n] = 0;
98  }
99 
100  int count = 0;
101  for (j = jmin; j <= jmax; j++) {
102  for (k = kmin; k <= kmax; k++) {
103  int itr = 0;
104  for (n = 0; n < body->nfacets; n++) {
105  if (cof[4][n] != 0) {
106  double yy, zz;
107  yy = body->surface[n].y3 - y[j];
108  zz = body->surface[n].z3 - z[k];
109  double l1, l2, l3;
110  l1 = (cof[1][n]*zz - cof[3][n]*yy) * cof[4][n];
111  l2 = (cof[2][n]*yy - cof[0][n]*zz) * cof[4][n];
112  l3 = 1 - l1 - l2;
113  if ((l1 > -eps) && (l2 > -eps) && (l3 > -eps)){
114  xd[itr] = l1*body->surface[n].x1 + l2*body->surface[n].x2 + l3*body->surface[n].x3;
115  dist[itr] = absolute(x[imin]-xd[itr]);
116  itr++;
117  }
118  }
119  }
120  if (itr > itr_max) {
121  if (!mpi->rank) {
122  fprintf(stderr,"Error: In IBIdentyBody() - itr > %d. Recompilation of code needed.\n",itr_max);
123  fprintf(stderr,"Increase the value of \"itr_max\" in IBIdentyBody().\n");
124  }
125  return(1);
126  }
127 
128  if (itr > 0) {
129  int ii, jj;
130 
131  for (ii = 0; ii < itr; ii++) {
132  for (jj = ii+1; jj < itr; jj++) {
133  if (dist[jj] < dist[ii]) {
134  double temp;
135  temp = dist[ii];
136  dist[ii] = dist[jj];
137  dist[jj] = temp;
138  temp = xd[ii];
139  xd[ii] = xd[jj];
140  xd[jj] = temp;
141  }
142  }
143  }
144 
145  for (ii = 1; ii < itr-1; ii++) {
146  if (absolute(xd[ii]-xd[ii-1]) < eps) {
147  for (jj = ii+1; jj < itr; jj++) {
148  xd[jj-1] = xd[jj];
149  dist[jj-1] = dist[jj];
150  }
151  itr--;
152  ii--;
153  }
154  }
155 
156  ii = 0;
157  int inside = 0;
158  for (i = imin; i <= imax; i++) {
159  if ((x[i]-xd[ii])*(x[i]-xd[ii+1]) < eps) {
160  inside = 1;
161  /* this point is inside */
162  /* check if (i,j,k) lies within this process */
163  /* if so, set blank */
164  if ( ((i-mpi->is[0])*(i-mpi->ie[0]) <= 0)
165  && ((j-mpi->is[1])*(j-mpi->ie[1]) <= 0)
166  && ((k-mpi->is[2])*(k-mpi->ie[2]) <= 0)) {
167  /* calculate local indices */
168  int index[_IB_NDIMS_];
169  index[0] = i-mpi->is[0];
170  index[1] = j-mpi->is[1];
171  index[2] = k-mpi->is[2];
172  int p; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,p);
173  blank[p] = 0;
174  count++;
175  }
176  } else {
177  if (inside) {
178  if (ii+2 < itr-1) ii += 2;
179  inside = 0;
180  }
181  }
182  }
183 
184  }
185  }
186  }
187 
188  free(xd);
189  free(dist);
190  return(count);
191 }
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
Body3D
Structure defining a body.
Definition: immersedboundaries.h:100
_IB_NDIMS_
#define _IB_NDIMS_
Definition: immersedboundaries.h:10
Body3D::surface
Facet3D * surface
Definition: immersedboundaries.h:102
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
absolute
#define absolute(a)
Definition: math_ops.h:32
max
#define max(a, b)
Definition: math_ops.h:18
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
min
#define min(a, b)
Definition: math_ops.h:14