immersedboundaries.h File Reference

Structures and function definitions for immersed boundaries. More...

#include <basic.h>

Go to the source code of this file.

Data Structures
struct	Facet3D
	Structure defining a facet. More...
struct	FacetMap
	Structure defining a facet map. More...
struct	Body3D
	Structure defining a body. More...
struct	IBNode
	Structure defining an immersed boundary node. More...
struct	ImmersedBoundary
	Structure containing variables for immersed boundary implementation. More...

Macros
#define	_IB_NDIMS_   3
#define	_IB_NNODES_   8
#define	_IB_XY_   "2d (xy)"
#define	_IB_XZ_   "2d (xz)"
#define	_IB_YZ_   "2d (yz)"
#define	_IB_3D_   "3d"

Functions
int	IBReadBodySTL (Body3D **, char *, void *, int *)
int	IBWriteBodySTL (Body3D *, char *, void *, int, int *)
int	IBCleanup (void *)
int	IBComputeBoundingBox (Body3D *)
int	IBCreateFacetMapping (void *, void *, double *, int *, int)
int	IBIdentifyBody (void *, int *, int *, int, void *, double *, double *)
int	IBIdentifyBoundary (void *, void *, int *, int, double *)
int	IBIdentifyMode (double *, int *, void *)
int	IBNearestFacetNormal (void *, void *, double *, double, int *, int)
int	IBInterpCoeffs (void *, void *, double *, int *, int, double *)
int	IBAssembleGlobalFacetData (void *, void *, const double *const, double **const, int)
int	IBComputeNormalGradient (void *, void *, const double *const, int, double **const)
int	IBComputeFacetVar (void *, void *, const double *const, int, double **const)

Detailed Description

Structures and function definitions for immersed boundaries.

Author

Debojyoti Ghosh

Definition in file immersedboundaries.h.

---

Data Structure Documentation

Facet3D

struct Facet3D

Structure defining a facet.

A "facet" is the basic unit of a tessellated surface in 3D: It is a triangle in an unstructured triangulated surface, defined by its vertices and surface normal.

See also

https://en.wikipedia.org/wiki/STL_(file_format)

Definition at line 39 of file immersedboundaries.h.

Data Fields
double	x1	x-coordinate of vertex 1
double	x2	x-coordinate of vertex 2
double	x3	x-coordinate of vertex 3
double	y1	y-coordinate of vertex 1
double	y2	y-coordinate of vertex 2
double	y3	y-coordinate of vertex 3
double	z1	z-coordinate of vertex 1
double	z2	z-coordinate of vertex 2
double	z3	z-coordinate of vertex 3
double	nx	x-component of surface normal
double	ny	y-component of surface normal
double	nz	z-component of surface normal

FacetMap

struct FacetMap

Structure defining a facet map.

A facet map contains information for facets that lie within the local computational domain of this MPI rank.

Definition at line 67 of file immersedboundaries.h.

Data Fields
Facet3D *	facet	pointer to the facet
int	index	index of this facet in the array Body3D::surface
int	interp_nodes[_IB_NNODES_]	indices of grid points surrounding the facet centroid
double	interp_coeffs[_IB_NNODES_]	interpolation coefficients corresponding to FacetMap::interp_nodes
int	interp_nodes_ns[_IB_NNODES_]	indices of grid points surrounding the "near-surface" point near the centroid
double	interp_coeffs_ns[_IB_NNODES_]	interpolation coefficients corresponding to FacetMap::interp_nodes_ns
double	xc	x-coordinate of centroid of FacetMap::facet
double	yc	y-coordinate of centroid of FacetMap::facet
double	zc	z-coordinate of centroid of FacetMap::facet
double	xns	x-coordinate of "near surface" point of FacetMap::facet
double	yns	y-coordinate of "near surface" point of FacetMap::facet
double	zns	z-coordinate of "near surface" point of FacetMap::facet
double	dx	FacetMap::xns - FacetMap::xc
double	dy	FacetMap::yns - FacetMap::yc
double	dz	FacetMap::zns - FacetMap::zc

Body3D

struct Body3D

Structure defining a body.

A 3D body whose surface is represented as a collection of faces of type Facet3D.

See also

https://en.wikipedia.org/wiki/STL_(file_format)

Definition at line 100 of file immersedboundaries.h.

Data Fields
int	nfacets	number of surface facets
Facet3D *	surface	array of surface facets
double	xmin	x-coordinate of lower end of bounding box
double	xmax	x-coordinate of higher end of bounding box
double	ymin	y-coordinate of lower end of bounding box
double	ymax	y-coordinate of higher end of bounding box
double	zmin	z-coordinate of lower end of bounding box
double	zmax	z-coordinate of higher end of bounding box

IBNode

struct IBNode

Structure defining an immersed boundary node.

An immersed boundary node is a grid point inside the immersed body but within stencil-width-distance of a point outside the body.

Definition at line 125 of file immersedboundaries.h.

Data Fields
int	i	grid index along x
int	j	grid index along y
int	k	grid index along z
int	p	array index of this point
double	x	x-coordinate of the boundary node
double	y	y-coordinate of the boundary node
double	z	z-coordinate of the boundary node
Facet3D *	face	the nearest facet of ImmersedBoundary::body
int	interp_nodes[_IB_NNODES_]	indices of interior nodes from which to extrapolate
double	interp_coeffs[_IB_NNODES_]	interpolation coefficients corresponding to IBNode::interp_nodes
double	interp_node_distance	Distance from the immersed body surface to the interior node from which to extrapolate
double	surface_distance	Distance from this node to the immersed body surface

ImmersedBoundary

struct ImmersedBoundary

Structure containing variables for immersed boundary implementation.

This structure contains all the variables needed to implement immersed boundaries.

Definition at line 153 of file immersedboundaries.h.

Data Fields
Body3D *	body	immersed body
IBNode *	boundary	immersed boundary nodes
FacetMap *	fmap	list of "local" facets
double	tolerance	zero tolerance
double	delta	small number
int	itr_max	maximum intersections in ray-tracing method
int	n_boundary_nodes	number of immersed boundary nodes
int	nfacets_local	number of "local" facets
char	mode[_MAX_STRING_SIZE_]	identifies if the simulation is 2D along a plane or truly 3D. See also IBIdentifyMode()

Macro Definition Documentation

_IB_NDIMS_

#define _IB_NDIMS_   3

Immersed boundaries are implemented only for 3-dimensional simulations.

Definition at line 10 of file immersedboundaries.h.

_IB_NNODES_

#define _IB_NNODES_   8

Number of grid points surrounding a random point in space, i.e., number of corners of a cube.

Definition at line 13 of file immersedboundaries.h.

_IB_XY_

#define _IB_XY_   "2d (xy)"

"Pseudo-2D" simulation in the x-y plane

Definition at line 16 of file immersedboundaries.h.

_IB_XZ_

#define _IB_XZ_   "2d (xz)"

"Pseudo-2D" simulation in the x-z plane

Definition at line 18 of file immersedboundaries.h.

_IB_YZ_

#define _IB_YZ_   "2d (yz)"

"Pseudo-2D" simulation in the y-z plane

Definition at line 20 of file immersedboundaries.h.

_IB_3D_

#define _IB_3D_   "3d"

3D simulation

Definition at line 22 of file immersedboundaries.h.

Function Documentation

IBReadBodySTL()

int IBReadBodySTL	(	Body3D **	body,
		char *	filename,
		void *	m,
		int *	stat
	)

Function to read a 3D surface from a STL file. See https://en.wikipedia.org/wiki/STL_(file_format) for details of the STL file format. Notes:

• This function only reads ASCII STL files.
• The body read from this file is distributed to all MPI ranks.

Parameters

body	3D body to be allocated and read from file
filename	Filename
m	MPI object of type MPIVariables
stat	Status (0: success; non-0: failure)

Definition at line 21 of file IBReadBodySTL.c.

{
  MPIVariables  *mpi = (MPIVariables*) m;
  int n, ierr;
  *stat = 0;

  if ((*body) != NULL) {
    if (!mpi->rank) {
      fprintf(stderr,"Error in IBReadBodySTL(): pointer to immersed body not NULL.\n");
    }
    *stat = -1;
    return(0);
  }

  /* Rank 0 reads in file */
  if (!mpi->rank) {

    FILE  *in;
    in = fopen(filename,"r");
    if (!in) {
      fprintf(stderr,"Error in IBReadBodySTL(): file %s not found or cannot be opened.\n",
              filename);
      *stat = 1;
    }
    else {

      int   nfacets = 0;
        char  word[_MAX_STRING_SIZE_];
      /* Count number of facets in STL file */
      ierr = fscanf(in,"%s",word);
      if (!strcmp(word, "solid")) {
        nfacets = 0;
        while (strcmp(word, "endsolid")) {
          ierr = fscanf(in,"%s",word);
          if (!strcmp(word, "facet")) nfacets++;
        }
      }
      fclose(in);

      if (nfacets > 0) {
        (*body) = (Body3D*) calloc (1,sizeof(Body3D));
        (*body)->nfacets = nfacets;
        (*body)->surface = (Facet3D*) calloc (nfacets,sizeof(Facet3D));

        /* Read in STL data */
        in = fopen(filename,"r");
        ierr = fscanf(in,"%s",word);
        n = 0;
        while (strcmp(word, "endsolid")) {
          double t1, t2, t3;
          ierr = fscanf(in,"%s",word);
          if (!strcmp(word, "facet")) {
            ierr = fscanf(in,"%s",word);
            if (strcmp(word, "normal")) {
              fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
            } else {
              ierr = fscanf(in,"%lf",&t1);
              ierr = fscanf(in,"%lf",&t2);
              ierr = fscanf(in,"%lf",&t3);
              (*body)->surface[n].nx = t1; 
              (*body)->surface[n].ny = t2; 
              (*body)->surface[n].nz = t3;
            }
            ierr = fscanf(in,"%s",word);
            if (strcmp(word, "outer")) {
              fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
              *stat = 1;
            } else {
              ierr = fscanf(in,"%s",word);
              if (strcmp(word, "loop")) {
                fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
                *stat = 1;
              } else {
                ierr = fscanf(in,"%s",word);
                if (strcmp(word, "vertex")) {
                  fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
                  *stat = 1;
                } else {
                  ierr = fscanf(in,"%lf",&t1);
                  ierr = fscanf(in,"%lf",&t2);
                  ierr = fscanf(in,"%lf",&t3);
                  (*body)->surface[n].x1 = t1; 
                  (*body)->surface[n].y1 = t2; 
                  (*body)->surface[n].z1 = t3;
                }
                ierr = fscanf(in,"%s",word);
                if (strcmp(word, "vertex")) {
                  fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
                  *stat = 1;
                } else {
                  ierr = fscanf(in,"%lf",&t1);
                  ierr = fscanf(in,"%lf",&t2);
                  ierr = fscanf(in,"%lf",&t3);
                  (*body)->surface[n].x2 = t1; 
                  (*body)->surface[n].y2 = t2; 
                  (*body)->surface[n].z2 = t3;
                }
                ierr = fscanf(in,"%s",word);
                if (strcmp(word, "vertex")) {
                  fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
                  *stat = 1;
                } else {
                  ierr = fscanf(in,"%lf",&t1);
                  ierr = fscanf(in,"%lf",&t2);
                  ierr = fscanf(in,"%lf",&t3);
                  (*body)->surface[n].x3 = t1; 
                  (*body)->surface[n].y3 = t2; 
                  (*body)->surface[n].z3 = t3;
                }
              }
            }
            n++;
          }
        }
        fclose(in);
        /* done reading STL file */
        if (n != nfacets) {
          fprintf(stderr,"Error in IBReadBodySTL(): Inconsistency in number of facets read.\n");
          free((*body)->surface);
          free(*body);
          *stat = 1;
        }
      } else {
        fprintf(stderr,"Error in IBReadBodySTL(): nfacets = 0!!\n");
        *stat = 1;
      }

    }
  }
  MPIBroadcast_integer(stat,1,0,&mpi->world);

  if ((*stat)) return(0);

  /* Distribute the body to all MPI ranks */
  int nfacets;
  if (!mpi->rank) nfacets = (*body)->nfacets;
  MPIBroadcast_integer(&nfacets,1,0,&mpi->world);

  if (mpi->rank) {
    (*body) = (Body3D*) calloc (1,sizeof(Body3D));
    (*body)->nfacets = nfacets;
    (*body)->surface = (Facet3D*) calloc (nfacets,sizeof(Facet3D));
  }

  int bufdim = 12;
  double *buffer = (double*) calloc (bufdim*nfacets,sizeof(double));
  if (!mpi->rank) {
    for (n=0; n<nfacets; n++) {
      buffer[n*bufdim+ 0] = (*body)->surface[n].x1;
      buffer[n*bufdim+ 1] = (*body)->surface[n].x2;
      buffer[n*bufdim+ 2] = (*body)->surface[n].x3;
      buffer[n*bufdim+ 3] = (*body)->surface[n].y1;
      buffer[n*bufdim+ 4] = (*body)->surface[n].y2;
      buffer[n*bufdim+ 5] = (*body)->surface[n].y3;
      buffer[n*bufdim+ 6] = (*body)->surface[n].z1;
      buffer[n*bufdim+ 7] = (*body)->surface[n].z2;
      buffer[n*bufdim+ 8] = (*body)->surface[n].z3;
      buffer[n*bufdim+ 9] = (*body)->surface[n].nx;
      buffer[n*bufdim+10] = (*body)->surface[n].ny;
      buffer[n*bufdim+11] = (*body)->surface[n].nz;
    }
  }
  MPIBroadcast_double(buffer,(nfacets*bufdim),0,&mpi->world);
  if (mpi->rank) {
    for (n=0; n<nfacets; n++) {
      (*body)->surface[n].x1 = buffer[n*bufdim+ 0];
      (*body)->surface[n].x2 = buffer[n*bufdim+ 1];
      (*body)->surface[n].x3 = buffer[n*bufdim+ 2];
      (*body)->surface[n].y1 = buffer[n*bufdim+ 3];
      (*body)->surface[n].y2 = buffer[n*bufdim+ 4];
      (*body)->surface[n].y3 = buffer[n*bufdim+ 5];
      (*body)->surface[n].z1 = buffer[n*bufdim+ 6];
      (*body)->surface[n].z2 = buffer[n*bufdim+ 7];
      (*body)->surface[n].z3 = buffer[n*bufdim+ 8];
      (*body)->surface[n].nx = buffer[n*bufdim+ 9];
      (*body)->surface[n].ny = buffer[n*bufdim+10];
      (*body)->surface[n].nz = buffer[n*bufdim+11];
    }
  }
  free(buffer);

  return(0);
}

IBWriteBodySTL()

int IBWriteBodySTL	(	Body3D *	body,
		char *	filename,
		void *	m,
		int	rank,
		int *	stat
	)

Function to write a 3D surface from a STL file. See https://en.wikipedia.org/wiki/STL_(file_format) for details of the STL file format. Notes:

• This function only writes ASCII STL files.

It is assumed that all MPI ranks have the body data. The MPI rank specified as the input rank will actually write the file.

Parameters

body	3D body to write to file
filename	Filename
m	MPI object of type MPIVariables
rank	MPI rank that does the writing
stat	Status (0: success; non-0: failure)

Definition at line 23 of file IBWriteBodySTL.c.

{
  MPIVariables  *mpi = (MPIVariables*) m;

  if (mpi->rank == rank) {
    FILE *out;
    out = fopen(filename,"w");
    if (!out) {
      fprintf(stderr,"Error in IBWriteBodySTL(): Unable to open file ");
      fprintf(stderr,"%s for writing on rank %d.\n",filename,mpi->rank);
      *stat = 1;
    } else {
      fprintf(out,"solid TEST\n");
      int n;
      for (n = 0; n < body->nfacets; n++) {
        fprintf(out,"  facet normal %1.16e %1.16e %1.16e\n",
                body->surface[n].nx,body->surface[n].ny,body->surface[n].nz);
        fprintf(out,"    outer loop\n");
        fprintf(out,"      vertex %1.16e %1.16e %1.16e\n",
                body->surface[n].x1,body->surface[n].y1,body->surface[n].z1);
        fprintf(out,"      vertex %1.16e %1.16e %1.16e\n",
                body->surface[n].x2,body->surface[n].y2,body->surface[n].z2);
        fprintf(out,"      vertex %1.16e %1.16e %1.16e\n",
                body->surface[n].x3,body->surface[n].y3,body->surface[n].z3);
        fprintf(out,"    endloop\n");
        fprintf(out,"  endfacet\n");
      }
      fprintf(out,"endsolid TEST\n");
      *stat = 0;
    }
  }
  return(0);
}

IBCleanup()

int IBCleanup

(

void *

)

Definition at line 11 of file IBCleanup.c.

{
  ImmersedBoundary *ib = (ImmersedBoundary*) s;
  if (!ib) return(0);

  free(ib->body->surface);
  free(ib->body);

  if (ib->n_boundary_nodes > 0) free(ib->boundary);
  if (ib->nfacets_local > 0) free(ib->fmap);

  return(0);
}

IBComputeBoundingBox()

int IBComputeBoundingBox

(

Body3D *

b

)

Compute the bounding box for a given body.

Parameters

b

The body

Definition at line 9 of file IBComputeBoundingBox.c.

{
  b->xmin = b->xmax = b->surface[0].x1;
  b->ymin = b->ymax = b->surface[0].y1;
  b->zmin = b->zmax = b->surface[0].z1;

  int n;
  for (n = 0; n < b->nfacets; n++) {
    if (b->surface[n].x1 < b->xmin) b->xmin = b->surface[n].x1;
    if (b->surface[n].x2 < b->xmin) b->xmin = b->surface[n].x2;
    if (b->surface[n].x3 < b->xmin) b->xmin = b->surface[n].x3;

    if (b->surface[n].y1 < b->ymin) b->ymin = b->surface[n].y1;
    if (b->surface[n].y2 < b->ymin) b->ymin = b->surface[n].y2;
    if (b->surface[n].y3 < b->ymin) b->ymin = b->surface[n].y3;

    if (b->surface[n].z1 < b->zmin) b->zmin = b->surface[n].z1;
    if (b->surface[n].z2 < b->zmin) b->zmin = b->surface[n].z2;
    if (b->surface[n].z3 < b->zmin) b->zmin = b->surface[n].z3;

    if (b->surface[n].x1 > b->xmax) b->xmax = b->surface[n].x1;
    if (b->surface[n].x2 > b->xmax) b->xmax = b->surface[n].x2;
    if (b->surface[n].x3 > b->xmax) b->xmax = b->surface[n].x3;

    if (b->surface[n].y1 > b->ymax) b->ymax = b->surface[n].y1;
    if (b->surface[n].y2 > b->ymax) b->ymax = b->surface[n].y2;
    if (b->surface[n].y3 > b->ymax) b->ymax = b->surface[n].y3;

    if (b->surface[n].z1 > b->zmax) b->zmax = b->surface[n].z1;
    if (b->surface[n].z2 > b->zmax) b->zmax = b->surface[n].z2;
    if (b->surface[n].z3 > b->zmax) b->zmax = b->surface[n].z3;
  }
  return(0);
}

IBCreateFacetMapping()

int IBCreateFacetMapping	(	void *	ib,
		void *	m,
		double *	X,
		int *	dim,
		int	ghosts
	)

This function creates a "facet map", i.e., on each MPI rank, it does the following:

• Makes a list of facets (defining the immersed body surface) that lie within the local computational domain of this MPI rank ("local facets").
• For each local facet, finds and stores the indices of the grid points that surround it, as well as the trilinear interpolation coefficients.
• For each local facet, finds a "near-surface" point, i.e., a point near the surface ("near" in terms of the local grid spacing) along the outward surface normal (i.e., outside the body), and finds and stores the indices of the grid points that surround it, as well as the trilinear interpolation coefficients.

Note: each MPI rank has a copy of the entire immersed body, i.e., all the facets.

Parameters

ib	Immersed boundary object of type ImmersedBoundary
m	MPI object of type MPIVariables
X	Array of local spatial coordinates
dim	Local dimensions
ghosts	Number of ghost points

Definition at line 143 of file IBCreateFacetMapping.c.

{
  ImmersedBoundary  *IB     = (ImmersedBoundary*) ib;
  MPIVariables      *mpi    = (MPIVariables*) m;
  Body3D            *body   = IB->body;
  int               nfacets = body->nfacets, n, count, ierr;
  Facet3D           *facets = body->surface;

  double  *x = X, 
          *y = (x + dim[0] + 2*ghosts), 
          *z = (y + dim[1] + 2*ghosts);

  double  xmin = 0.5 * (x[ghosts-1]         + x[ghosts]),
          xmax = 0.5 * (x[dim[0]+ghosts-1]  + x[dim[0]+ghosts]),
          ymin = 0.5 * (y[ghosts-1]         + y[ghosts]),
          ymax = 0.5 * (y[dim[1]+ghosts-1]  + y[dim[1]+ghosts]),
          zmin = 0.5 * (z[ghosts-1]         + z[ghosts]),
          zmax = 0.5 * (z[dim[2]+ghosts-1]  + z[dim[2]+ghosts]);

  count = 0;
  for (n = 0; n < nfacets; n++) {

    /* find facet centroid */
    double xc, yc, zc;
    xc = (facets[n].x1 + facets[n].x2 + facets[n].x3) / 3.0;
    yc = (facets[n].y1 + facets[n].y2 + facets[n].y3) / 3.0;
    zc = (facets[n].z1 + facets[n].z2 + facets[n].z3) / 3.0;

    if (!strcmp(IB->mode,_IB_3D_)) {
      if (isInside(xc,xmin,xmax) && isInside(yc,ymin,ymax) && isInside(zc,zmin,zmax)) count++;
    } else if (!strcmp(IB->mode,_IB_XY_)) {
      if (isInside(xc,xmin,xmax) && isInside(yc,ymin,ymax)) count++;
    } else if (!strcmp(IB->mode,_IB_XZ_)) {
      if (isInside(xc,xmin,xmax) && isInside(zc,zmin,zmax)) count++;
    } else if (!strcmp(IB->mode,_IB_YZ_)) {
      if (isInside(yc,ymin,ymax) && isInside(zc,zmin,zmax)) count++;
    }
  }

  int nfacets_local = count;
  if (nfacets_local > 0) {

    FacetMap *fmap = (FacetMap*) calloc (nfacets_local, sizeof(FacetMap));
    count = 0;

    for (n = 0; n < nfacets; n++) {

      double xc, yc, zc;
      xc = (facets[n].x1 + facets[n].x2 + facets[n].x3) / 3.0;
      yc = (facets[n].y1 + facets[n].y2 + facets[n].y3) / 3.0;
      zc = (facets[n].z1 + facets[n].z2 + facets[n].z3) / 3.0;

      int flag = 0;
      if (!strcmp(IB->mode,_IB_3D_)) {
        if (isInside(xc,xmin,xmax) && isInside(yc,ymin,ymax) && isInside(zc,zmin,zmax)) flag = 1;
      } else if (!strcmp(IB->mode,_IB_XY_)) {
        if (isInside(xc,xmin,xmax) && isInside(yc,ymin,ymax)) flag = 1;
      } else if (!strcmp(IB->mode,_IB_XZ_)) {
        if (isInside(xc,xmin,xmax) && isInside(zc,zmin,zmax)) flag = 1;
      } else if (!strcmp(IB->mode,_IB_YZ_)) {
        if (isInside(yc,ymin,ymax) && isInside(zc,zmin,zmax)) flag = 1;
      }

      if (flag == 1) {
        fmap[count].facet = facets + n;
        fmap[count].index = n;
        
        fmap[count].xc = xc;
        fmap[count].yc = yc;
        fmap[count].zc = zc;

        int ic,jc, kc;

        ierr = interpNodesCoeffs( mpi,
                                  xc, yc, zc,
                                  x, y, z,
                                  dim,
                                  ghosts,
                                  IB->mode,
                                  &ic, &jc, &kc,
                                  fmap[count].interp_nodes,
                                  fmap[count].interp_coeffs );
        if (ierr) {
          fprintf(stderr, "Error in IBCreateFacetMapping(): \n");
          fprintf(stderr, "  interpNodesCoeffs() returned with error code %d on rank %d.\n",
                  ierr, mpi->rank );
          return(ierr);
        }

        double dx = x[ic] - x[ic-1];
        double dy = y[jc] - y[jc-1];
        double dz = z[kc] - z[kc-1];
        double ds;
        if      (!strcmp(IB->mode,_IB_XY_)) ds = min(dx,dy);
        else if (!strcmp(IB->mode,_IB_XZ_)) ds = min(dx,dz);
        else if (!strcmp(IB->mode,_IB_YZ_)) ds = min(dy,dz);
        else                                ds = min3(dx,dy,dz);

        double nx = fmap[count].facet->nx;
        double ny = fmap[count].facet->ny;
        double nz = fmap[count].facet->nz;

        if (nx == 0.0) nx += IB->delta*ds;
        if (ny == 0.0) ny += IB->delta*ds;
        if (nz == 0.0) nz += IB->delta*ds;

        double xns = xc + sign(nx)*ds;
        double yns = yc + sign(ny)*ds;
        double zns = zc + sign(nz)*ds;

        fmap[count].dx = xns - xc;
        fmap[count].dy = yns - yc;
        fmap[count].dz = zns - zc;

        ierr = interpNodesCoeffs( mpi,
                                  xns, yns, zns,
                                  x, y, z,
                                  dim,
                                  ghosts,
                                  IB->mode,
                                  NULL,NULL,NULL,
                                  fmap[count].interp_nodes_ns,
                                  fmap[count].interp_coeffs_ns );
        if (ierr) {
          fprintf(stderr, "Error in IBCreateFacetMapping(): \n");
          fprintf(stderr, "  interpNodesCoeffs() returned with error code %d on rank %d.\n",
                  ierr, mpi->rank );
          return(ierr);
        }

        count++;
      }
    }

    IB->nfacets_local = nfacets_local;
    IB->fmap = fmap;  

  } else {

    IB->nfacets_local = 0;
    IB->fmap = NULL;

  }

  return(0);
}

IBIdentifyBody()

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

ib	Immersed boundary object of type ImmersedBoundary
dim_g	global dimensions
dim_l	local dimensions
ghosts	number of ghost points
m	MPI object of type MPIVariables
X	Array of global spatial coordinates
blank	Blanking array: for grid points within the body, this value will be set to 0

Definition at line 20 of file IBIdentifyBody.c.

{
  ImmersedBoundary  *IB     = (ImmersedBoundary*) ib;
  MPIVariables      *mpi    = (MPIVariables*) m;
  Body3D            *body   = IB->body;
  double            *x      = X, 
                    *y      = X+dim_g[0], 
                    *z      = X+dim_g[0]+dim_g[1],
                    eps     = IB->tolerance;
  int               itr_max = IB->itr_max,
                      i, j, k, n, v;

    double xmax = body->xmax, 
         xmin = body->xmin, 
         ymax = body->ymax, 
         ymin = body->ymin, 
         zmax = body->zmax, 
         zmin = body->zmin;
  
  double fac = 1.5;
    double Lx, Ly, Lz;
    Lx = (xmax - xmin);
    Ly = (ymax - ymin);
    Lz = (zmax - zmin);
    double xc, yc, zc;
    xc = 0.5 * (xmin + xmax);
    yc = 0.5 * (ymin + ymax);
    zc = 0.5 * (zmin + zmax);
    xmax = xc + fac * Lx/2;
    xmin = xc - fac * Lx/2;
    ymax = yc + fac * Ly/2;
    ymin = yc - fac * Ly/2;
    zmax = zc + fac * Lz/2;
    zmin = zc - fac * Lz/2;

    int imin, imax, jmin, jmax, kmin, kmax;
    imin = dim_g[0]-1;  imax = 0;
    jmin = dim_g[1]-1;  jmax = 0;
    kmin = dim_g[2]-1;  kmax = 0;
    for (i = 0; i < dim_g[0]; i++) {
        for (j = 0; j < dim_g[1]; j++) {
            for (k = 0; k < dim_g[2]; k++) {
                if (   ((x[i]-xmin)*(x[i]-xmax) < 0) 
            && ((y[j]-ymin)*(y[j]-ymax) < 0) 
            && ((z[k]-zmin)*(z[k]-zmax) < 0)) {
                    imin = min(i, imin);
                    imax = max(i, imax);
                    jmin = min(j, jmin);
                    jmax = max(j, jmax);
                    kmin = min(k, kmin);
                    kmax = max(k, kmax);
                }
            }
        }
    }

    double *cof[5], *xd, *dist;
  xd    = (double*) calloc (itr_max,sizeof(double));
  dist  = (double*) calloc (itr_max,sizeof(double));
    for (v = 0; v < 5; v++) cof[v] = (double*) calloc(body->nfacets, sizeof(double));

    for (n = 0; n < body->nfacets; n++) {
        cof[0][n] = body->surface[n].y1 - body->surface[n].y3;
        cof[1][n] = body->surface[n].y2 - body->surface[n].y3;
    cof[2][n] = body->surface[n].z1 - body->surface[n].z3;
    cof[3][n] = body->surface[n].z2 - body->surface[n].z3;
        double den = cof[0][n]*cof[3][n] - cof[1][n]*cof[2][n];
        if (absolute(den) > eps)    cof[4][n] = 1.0 / den;
        else                            cof[4][n] = 0;      
    }

  int count = 0;
    for (j = jmin; j <= jmax; j++) {
        for (k = kmin; k <= kmax; k++) {
            int itr = 0;
            for (n = 0; n < body->nfacets; n++) {
                if (cof[4][n] != 0) {
                    double yy, zz;
                    yy = body->surface[n].y3 - y[j];
                    zz = body->surface[n].z3 - z[k];
                    double l1, l2, l3;
                    l1 = (cof[1][n]*zz - cof[3][n]*yy) * cof[4][n];
                    l2 = (cof[2][n]*yy - cof[0][n]*zz) * cof[4][n];
                    l3 = 1 - l1 - l2;
                    if ((l1 > -eps) && (l2 > -eps) && (l3 > -eps)){
                        xd[itr]   = l1*body->surface[n].x1 + l2*body->surface[n].x2 + l3*body->surface[n].x3;
                        dist[itr] = absolute(x[imin]-xd[itr]);
                        itr++;
                    }
                }
            }
            if (itr > itr_max) {
        if (!mpi->rank) {
                  fprintf(stderr,"Error: In IBIdentyBody() - itr > %d. Recompilation of code needed.\n",itr_max);
                  fprintf(stderr,"Increase the value of \"itr_max\" in IBIdentyBody().\n");
        }
        return(1);
            }
    
            if (itr > 0) {
                int ii, jj;
    
                for (ii = 0; ii < itr; ii++) {
                    for (jj = ii+1; jj < itr; jj++) {
                        if (dist[jj] < dist[ii]) {
                            double temp;
                            temp     = dist[ii];
                            dist[ii] = dist[jj];
                            dist[jj] = temp;
                            temp     = xd[ii];
                            xd[ii]   = xd[jj];
                            xd[jj]   = temp;
                        }
                    }
                }

                for (ii = 1; ii < itr-1; ii++) {
                    if (absolute(xd[ii]-xd[ii-1]) < eps) {
                        for (jj = ii+1; jj < itr; jj++) {
                            xd[jj-1]    = xd[jj];
                            dist[jj-1]  = dist[jj];
                        }
                        itr--;
                        ii--;
                    }
                }

                ii = 0;
                int inside = 0;
                for (i = imin; i <= imax; i++) {
                    if ((x[i]-xd[ii])*(x[i]-xd[ii+1]) < eps) {
                        inside = 1;
            /* this point is inside                      */
            /* check if (i,j,k) lies within this process */
            /* if so, set blank                          */
            if (   ((i-mpi->is[0])*(i-mpi->ie[0]) <= 0) 
                && ((j-mpi->is[1])*(j-mpi->ie[1]) <= 0) 
                && ((k-mpi->is[2])*(k-mpi->ie[2]) <= 0)) {
              /* calculate local indices */
              int index[_IB_NDIMS_];
              index[0] = i-mpi->is[0];
              index[1] = j-mpi->is[1];
              index[2] = k-mpi->is[2];
              int p; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,p);
              blank[p] = 0;
              count++;
            }
                    } else {
                        if (inside) {
                            if (ii+2 < itr-1)   ii += 2;
                            inside = 0;
                        }
                    }
                }

            }
        }
    }

  free(xd);
  free(dist);
  return(count);
}

IBIdentifyBoundary()

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

Identify the immersed boundary points: an immersed boundary point is any grid point inside the immersed body that is within stencil-width-distance of a grid point outside the immersed body. This function does the following:

• count the number of immersed boundary points.
• allocate the array of immersed boundary points and set their indices.

Parameters

ib	Immersed boundary object of type ImmersedBoundary
m	MPI object of type MPIVariables
dim_l	local dimensions
ghosts	number of ghost points
blank	Blanking array: for grid points within the immersed body, this value will be set to 0

Definition at line 158 of file IBIdentifyBoundary.c.

{
  ImmersedBoundary  *IB     = (ImmersedBoundary*) ib;
  MPIVariables      *mpi    = (MPIVariables*) m;
  Body3D            *body   = IB->body;

  int imax = dim_l[0],
      jmax = dim_l[1],
      kmax = dim_l[2];

  int n_boundary_nodes = CountBoundaryPoints(imax,jmax,kmax,ghosts,blank);
  IB->n_boundary_nodes = n_boundary_nodes;
  if (n_boundary_nodes == 0) IB->boundary = NULL;
  else {
    IB->boundary = (IBNode*) calloc (n_boundary_nodes, sizeof(IBNode));
    int check = SetBoundaryPoints(imax,jmax,kmax,ghosts,blank,IB->boundary);
    if (check != n_boundary_nodes) {
      fprintf(stderr,"Error in IBIdentifyBoundary(): Inconsistency encountered when setting boundary indices. ");
      fprintf(stderr,"on rank %d.\n",mpi->rank);
      fprintf(stderr,"SetBoundaryPoints() returned %d, while n_boundary_nodes is %d.\n",check,n_boundary_nodes);
    }
  }

  return(n_boundary_nodes);
}

IBIdentifyMode()

int IBIdentifyMode	(	double *	X,
		int *	dim,
		void *	ib
	)

Identify the "mode", i.e., whether the simulation is a true 3D simulation, or a 2D simulation being run in 3D. If extent of the immersed body is larger than the grid along a particular axis (say, x), then we assume that the intention is to simulate a 2D case around a 2D body in the plane normal to that axis ( y-z plane ).

For example, to simulate a 2D cylinder in the x-y plane, we consider a cylinder whose extent along z is larger than the extent of the grid along z, i.e., it sticks out of the computational domain at both ends.

If the immersed body is completely contained within the computational domain, or sticks out only on one end along a particular dimension, we assume it's a 3D simulation.

Parameters

X	Array of global spatial coordinates
dim	global dimensions
ib	Immersed boundary object of type ImmersedBoundary

Definition at line 24 of file IBIdentifyMode.c.

{
  ImmersedBoundary  *IB = (ImmersedBoundary*) ib;
  Body3D            *body   = IB->body;

  double *x = X, *y = x + dim[0], *z = y + dim[1];

  double grid_xmin = x[0],
         grid_xmax = x[dim[0]-1],
         grid_ymin = y[0],
         grid_ymax = y[dim[1]-1],
         grid_zmin = z[0],
         grid_zmax = z[dim[2]-1];

  if      ( (grid_xmin > body->xmin) && (grid_xmax < body->xmax) ) strcpy(IB->mode,_IB_YZ_);
  else if ( (grid_ymin > body->ymin) && (grid_ymax < body->ymax) ) strcpy(IB->mode,_IB_XZ_);
  else if ( (grid_zmin > body->zmin) && (grid_zmax < body->zmax) ) strcpy(IB->mode,_IB_XY_);
  else                                                             strcpy(IB->mode,_IB_3D_);

  return(0);
}

IBNearestFacetNormal()

int IBNearestFacetNormal	(	void *	ib,
		void *	m,
		double *	X,
		double	large_distance,
		int *	dim_l,
		int	ghosts
	)

For each immersed boundary point, find the nearest facet (Facet3D) of the immersed body (ImmersedBoundary::body). The "nearest" facet is the one which is closest to the boundary point in terms of the distance along the normal defined for that facet.

• The function will first try to find the nearest facet for which a line starting from the boundary point along the direction defind by the facet normal passes through that facet.
• Failing the above criterion, the function will find the nearest facet.

Note: This function is sensitive to the fact that the normals defined for the immersed body are the outward normals, i.e., pointing away from the body. If this function returns an error, make sure this is true in the STL file the body is defined in.

Parameters

ib	Immersed boundary object of type ImmersedBoundary
m	MPI object of type MPIVariables
X	Array of (local) spatial coordinates
large_distance	A large distance
dim_l	Integer array of local grid size in each spatial dimension
ghosts	Number of ghost points

Definition at line 25 of file IBNearestFacetNormal.c.

{
  ImmersedBoundary  *IB       = (ImmersedBoundary*) ib;
  MPIVariables      *mpi      = (MPIVariables*) m;
  Body3D            *body     = IB->body;
  Facet3D           *surface  = body->surface;
  IBNode            *boundary = IB->boundary;

  double  eps = IB->tolerance;
  int     nb  = IB->n_boundary_nodes,
          nf  = body->nfacets;

    int i, j, k, dg, n;
    for (dg = 0; dg < nb; dg++) {
        i = boundary[dg].i;
        j = boundary[dg].j;
        k = boundary[dg].k;

    double xp, yp, zp;
    _GetCoordinate_(0,i,dim_l,ghosts,X,xp);
    _GetCoordinate_(1,j,dim_l,ghosts,X,yp);
    _GetCoordinate_(2,k,dim_l,ghosts,X,zp);
    boundary[dg].x = xp;
    boundary[dg].y = yp;
    boundary[dg].z = zp;

        double  dist_min = large_distance;
        int     n_min    = -1;

        for (n = 0; n < nf; n++) {
            double x1, x2, x3;
            double y1, y2, y3;
            double z1, z2, z3;
            x1 = surface[n].x1; x2 = surface[n].x2; x3 = surface[n].x3;
            y1 = surface[n].y1; y2 = surface[n].y2; y3 = surface[n].y3;
            z1 = surface[n].z1; z2 = surface[n].z2; z3 = surface[n].z3;
            double dist =   surface[n].nx*(xp-surface[n].x1) 
                    + surface[n].ny*(yp-surface[n].y1) 
                    + surface[n].nz*(zp-surface[n].z1);
            if (dist > 0)   continue;
            if (absolute(dist) < dist_min) {
                short   is_it_in = 0;
                double  x_int, y_int, z_int;
                x_int = xp - dist * surface[n].nx;
                y_int = yp - dist * surface[n].ny;
                z_int = zp - dist * surface[n].nz;
                if (absolute(surface[n].nx) > eps) {
                    double den = (z2-z3)*(y1-y3)-(y2-y3)*(z1-z3);
                    double l1, l2, l3;
                    l1 = ((y2-y3)*(z3-z_int)-(z2-z3)*(y3-y_int)) / den;
                    l2 = ((z1-z3)*(y3-y_int)-(y1-y3)*(z3-z_int)) / den;
                    l3 = 1 - l1 - l2;
                    if ((l1 > -eps) && (l2 > -eps) && (l3 > -eps))  is_it_in = 1;
                } else if (absolute(surface[n].ny) > eps) {
                    double den = (x2-x3)*(z1-z3)-(z2-z3)*(x1-x3);
                    double l1, l2, l3;
                    l1 = ((z2-z3)*(x3-x_int)-(x2-x3)*(z3-z_int)) / den;
                    l2 = ((x1-x3)*(z3-z_int)-(z1-z3)*(x3-x_int)) / den;
                    l3 = 1 - l1 - l2;
                    if ((l1 > -eps) && (l2 > -eps) && (l3 > -eps))  is_it_in = 1;
                } else {
                    double den = (y2-y3)*(x1-x3)-(x2-x3)*(y1-y3);
                    double l1, l2, l3;
                    l1 = ((x2-x3)*(y3-y_int)-(y2-y3)*(x3-x_int)) / den;
                    l2 = ((y1-y3)*(x3-x_int)-(x1-x3)*(y3-y_int)) / den;
                    l3 = 1 - l1 - l2;
                    if ((l1 > -eps) && (l2 > -eps) && (l3 > -eps))  is_it_in = 1;
                }
                if (is_it_in) {
                    dist_min = absolute(dist);
                    n_min = n; 
                }
            }
        }
        if (n_min == -1) {
            for (n = 0; n < nf; n++) {
                double dist =   surface[n].nx*(xp-surface[n].x1) 
                      + surface[n].ny*(yp-surface[n].y1) 
                      + surface[n].nz*(zp-surface[n].z1);
                if (dist > eps) continue;
                else {
                    if (absolute(dist) < dist_min) {
                        dist_min = absolute(dist);
                        n_min = n;
                    }
                }
            }
        }

        if (n_min == -1)    {
      fprintf(stderr,"Error in IBNearestFacetNormal(): no nearest normal found for boundary node (%d,%d,%d) ",i,j,k);
      fprintf(stderr,"on rank %d.\n",mpi->rank);
      return(1);
    } else boundary[dg].face = &surface[n_min];
    }

  return(0);
}

IBInterpCoeffs()

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

ib	Immersed boundary object of type ImmersedBoundary
m	MPI object of type MPIVariables
X	Array of (local) spatial coordinates
dim_l	Integer array of local grid size in each spatial dimension
ghosts	Number of ghost points
blank	Blanking array: for grid points within the body, this value will be set to 0

Definition at line 22 of file IBInterpCoeffs.c.

{
  ImmersedBoundary  *IB       = (ImmersedBoundary*) ib;
  MPIVariables      *mpi      = (MPIVariables*) m;
  IBNode            *boundary = IB->boundary;

  double  eps         = IB->tolerance;
  int     maxiter     = IB->itr_max,
          n_boundary  = IB->n_boundary_nodes;

  int imax        = dim_l[0],
      jmax        = dim_l[1],
      kmax        = dim_l[2];

  int is = mpi->is[0],
      js = mpi->is[1],
      ks = mpi->is[2];

  static int index[_IB_NDIMS_];

    int dg;
    for (dg = 0; dg < n_boundary; dg++) {
        int    i, j, k, p;
        double xb, yb, zb;
        double nx, ny, nz;
        double xx, yy, zz;
        double dx, dy, dz;
    double ds, dist;
        double xtip, ytip, ztip;

        index[0] = i = boundary[dg].i;
        index[1] = j = boundary[dg].j;
        index[2] = k = boundary[dg].k;
    _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,p);

        xb = boundary[dg].x;
        yb = boundary[dg].y;
        zb = boundary[dg].z;

        nx = boundary[dg].face->nx;
        ny = boundary[dg].face->ny;
        nz = boundary[dg].face->nz;
        xx = boundary[dg].face->x1;
        yy = boundary[dg].face->y1;
        zz = boundary[dg].face->z1;

        dist = nx*(xx-xb) + ny*(yy-yb) + nz*(zz-zb);

    double x1, x2, y1, y2, z1, z2;
    _GetCoordinate_(0,(i+1),dim_l,ghosts,X,x1);
    _GetCoordinate_(0,(i-1),dim_l,ghosts,X,x2);
    _GetCoordinate_(1,(j+1),dim_l,ghosts,X,y1);
    _GetCoordinate_(1,(j-1),dim_l,ghosts,X,y2);
    _GetCoordinate_(2,(k+1),dim_l,ghosts,X,z1);
    _GetCoordinate_(2,(k-1),dim_l,ghosts,X,z2);
    dx = 0.5 * (x1 - x2);
    dy = 0.5 * (y1 - y2);
    dz = 0.5 * (z1 - z2);
        ds = min3(dx, dy, dz);

        xtip = xb + dist*nx;
        ytip = yb + dist*ny;
        ztip = zb + dist*nz;

        int is_it_in = 0;
        int iter = 0;
        int itip, jtip, ktip;
        while(!is_it_in && (iter < maxiter)) {
            iter++;
            itip = i;
            jtip = j;
            ktip = k;
          
            if (xtip > xb)  {
        double xx;
        _GetCoordinate_(0,itip,dim_l,ghosts,X,xx);
        while ((xx < xtip) && (itip < imax+ghosts-1)) {
          itip++;
          _GetCoordinate_(0,itip,dim_l,ghosts,X,xx);
        }
      } else {
        double xx;
        _GetCoordinate_(0,(itip-1),dim_l,ghosts,X,xx);
        while ((xx > xtip) && (itip > -ghosts)) {
          itip--;
          _GetCoordinate_(0,(itip-1),dim_l,ghosts,X,xx);
        }
      }

            if (ytip > yb) {
        double yy;
        _GetCoordinate_(1,jtip,dim_l,ghosts,X,yy);
        while ((yy < ytip) && (jtip < jmax+ghosts-1)) {
          jtip++;
          _GetCoordinate_(1,jtip,dim_l,ghosts,X,yy);
        }
      } else {
        double yy;
        _GetCoordinate_(1,(jtip-1),dim_l,ghosts,X,yy);
        while ((yy > ytip) && (jtip > -ghosts)) {
          jtip--;
          _GetCoordinate_(1,(jtip-1),dim_l,ghosts,X,yy);
        }
      }

            if (ztip > zb) {
        double zz;
        _GetCoordinate_(2,ktip,dim_l,ghosts,X,zz);
        while ((zz < ztip) && (ktip < kmax+ghosts-1))   {
          ktip++;
          _GetCoordinate_(2,ktip,dim_l,ghosts,X,zz);
        }
      } else {
        double zz;
        _GetCoordinate_(2,(ktip-1),dim_l,ghosts,X,zz);
        while ((zz > ztip) && (ktip > -ghosts)) {
          ktip--;
          _GetCoordinate_(2,(ktip-1),dim_l,ghosts,X,zz);
        }
      }

      int ptip[_IB_NNODES_];
      index[0] = itip  ; index[1] = jtip  ; index[2] = ktip  ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[0]);
      index[0] = itip-1; index[1] = jtip  ; index[2] = ktip  ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[1]);
      index[0] = itip  ; index[1] = jtip-1; index[2] = ktip  ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[2]);
      index[0] = itip  ; index[1] = jtip  ; index[2] = ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[3]);
      index[0] = itip-1; index[1] = jtip-1; index[2] = ktip  ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[4]);
      index[0] = itip  ; index[1] = jtip-1; index[2] = ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[5]);
      index[0] = itip-1; index[1] = jtip  ; index[2] = ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[6]);
      index[0] = itip-1; index[1] = jtip-1; index[2] = ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[7]);

            int nflow = 0;
            nflow += blank[ptip[0]];
            nflow += blank[ptip[1]];
            nflow += blank[ptip[2]];
            nflow += blank[ptip[3]];
            nflow += blank[ptip[4]];
            nflow += blank[ptip[5]];
            nflow += blank[ptip[6]];
            nflow += blank[ptip[7]];
            if (nflow == _IB_NNODES_) {
                is_it_in = 1;
            } else if (nflow < _IB_NNODES_) {
                is_it_in = 0;
                xtip += nx*absolute(ds);
                ytip += ny*absolute(ds);
                ztip += nz*absolute(ds);
            } else {
                fprintf(stderr,"Error in IBInterpCoeffs() (Bug in code) - counting interior points surrounding probe tip \n");
        fprintf(stderr,"on rank %d.\n", mpi->rank);
        fprintf(stderr,"Value of nflow is %d but can only be positive and <= %d.\n",nflow,_IB_NNODES_);
        return(1);
            }
        }

    if (!is_it_in) {
      fprintf(stderr,"Error in IBInterpCoeffs() on rank %d - interior point not found for immersed boundary point (%d,%d,%d)!\n",
              mpi->rank, i, j, k);
      return(1);
    }    
    
        double tlx[2],tly[2],tlz[2];
    _GetCoordinate_(0,(itip-1),dim_l,ghosts,X,tlx[0]);
    _GetCoordinate_(0,(itip  ),dim_l,ghosts,X,tlx[1]);
    _GetCoordinate_(1,(jtip-1),dim_l,ghosts,X,tly[0]);
    _GetCoordinate_(1,(jtip  ),dim_l,ghosts,X,tly[1]);
    _GetCoordinate_(2,(ktip-1),dim_l,ghosts,X,tlz[0]);
    _GetCoordinate_(2,(ktip  ),dim_l,ghosts,X,tlz[1]);

    int ptip[_IB_NNODES_];
    index[0]=itip-1; index[1]=jtip-1; index[2]=ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[0]);
    index[0]=itip  ; index[1]=jtip-1; index[2]=ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[1]);
    index[0]=itip-1; index[1]=jtip  ; index[2]=ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[2]);
    index[0]=itip  ; index[1]=jtip  ; index[2]=ktip-1; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[3]);
    index[0]=itip-1; index[1]=jtip-1; index[2]=ktip  ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[4]);
    index[0]=itip  ; index[1]=jtip-1; index[2]=ktip  ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[5]);
    index[0]=itip-1; index[1]=jtip  ; index[2]=ktip  ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[6]);
    index[0]=itip  ; index[1]=jtip  ; index[2]=ktip  ; _ArrayIndex1D_(_IB_NDIMS_,dim_l,index,ghosts,ptip[7]);
    _ArrayCopy1D_(ptip,boundary[dg].interp_nodes,_IB_NNODES_);

    double coeffs[_IB_NNODES_];
    TrilinearInterpCoeffs(tlx[0],tlx[1],tly[0],tly[1],tlz[0],tlz[1],xtip,ytip,ztip,&coeffs[0]);
    _ArrayCopy1D_(coeffs,boundary[dg].interp_coeffs,_IB_NNODES_);

        double tipdist = absolute(nx*(xx-xtip) + ny*(yy-ytip) + nz*(zz-ztip));
    boundary[dg].interp_node_distance = tipdist;
    boundary[dg].surface_distance = absolute(dist);
        if (tipdist < eps) {
            fprintf(stderr,"Warning in IBInterpCoeffs() on rank %d - how can probe tip be on surface? Tipdist = %e\n",
              mpi->rank,tipdist);
        }

    }
  return(0);
}

IBAssembleGlobalFacetData()

int IBAssembleGlobalFacetData	(	void *	m,
		void *	ib,
		const double *const	local_var,
		double **const	global_var,
		int	nvars
	)

Assemble any local data computed at immersed body surface (facet centroids) into a global array.

The local array must be of length (ImmersedBoundary::nfacets_local X nvars).

The global array must be NULL at input; it will have the size (Body3D::nfacets X nvars) at output on rank 0; it will remain NULL on other ranks.

Parameters

m	MPI object of type MPIVariables
ib	Immersed boundary object of type ImmersedBoundary
local_var	Local array
global_var	Array to store the gradient; must be NULL at input
nvars	Number of components in var

Definition at line 22 of file IBAssembleGlobalFacetData.c.

{
  MPIVariables     *mpi = (MPIVariables*) m; 
  ImmersedBoundary *IB  = (ImmersedBoundary*) ib;

  if ((*global_var) != NULL) {
    fprintf(stderr,"Error in IBAssembleGlobalFacetData()\n");
    fprintf(stderr," global_var is not NULL on rank %d\n",
            mpi->rank );
    return 1;
  }

  int nfacets_local = IB->nfacets_local;
  FacetMap *fmap = IB->fmap;

  if ((nfacets_local == 0) && (local_var != NULL)) {
    fprintf(stderr, "Error in IBAssembleGlobalFacetData()\n");
    fprintf(stderr, " nfacets_local is 0 but local_var is not NULL!\n");
    return 1;
  }
  if ((nfacets_local > 0) && (local_var == NULL)) {
    fprintf(stderr, "Error in IBAssembleGlobalFacetData()\n");
    fprintf(stderr, " nfacets_local > 0 but local_var is NULL!\n");
    return 1;
  }

#ifndef serial
  MPI_Status status;
#endif

    if (!mpi->rank) {

    int nfacets_global = IB->body->nfacets;

    /* allocate arrays for whole domain */
        *global_var = (double*) calloc (nfacets_global*nvars, sizeof(double));
    _ArraySetValue_((*global_var), (nfacets_global*nvars), 0.0);

    /* check array - to avoid double counting of facets */
        int *check = (int*) calloc (nfacets_global, sizeof(int));
    _ArraySetValue_(check, nfacets_global, 0);

    /* local data */
        for (int n = 0; n < nfacets_local; n++) {
      _ArrayAXPY_((local_var+n), 1.0, ((*global_var)+fmap[n].index), nvars);
          check[fmap[n].index]++;
        }

#ifndef serial
        for (int proc = 1; proc < mpi->nproc; proc++) {

            int nf_incoming;
            MPI_Recv(&nf_incoming, 1, MPI_INT, proc, 98927, MPI_COMM_WORLD, &status);

      if (nf_incoming > 0) {

              int    *indices_incoming = (int*)    calloc(nf_incoming, sizeof(int));
              double *var_incoming     = (double*) calloc(nf_incoming*nvars, sizeof(double));

              MPI_Recv(indices_incoming, nf_incoming, MPI_INT, proc, 98928, mpi->world, &status);
              MPI_Recv(var_incoming, nf_incoming*nvars, MPI_DOUBLE, proc, 98929, mpi->world, &status);

              for (int n = 0; n < nf_incoming; n++) {
          _ArrayAXPY_((var_incoming+n), 1.0, ((*global_var)+indices_incoming[n]), nvars);
                  check[indices_incoming[n]]++;
              }
        
        free(indices_incoming);
        free(var_incoming);
      }
        }
#endif

    for (int n = 0; n < nfacets_global; n++) {
      if (check[n] == 0) {
        fprintf(stderr, "Error in IBAssembleGlobalFacetData()\n");
        fprintf(stderr, "  No data received for facet %d\n", n);
        return 1;
      } else {
        _ArrayScale1D_(((*global_var)+n), (1.0/(double)check[n]), nvars);
      }
    }

    } else {

#ifndef serial

        MPI_Send(&nfacets_local, 1, MPI_INT, 0, 98927, MPI_COMM_WORLD);

    if (nfacets_local > 0) {

      int i, *indices = (int*) calloc (nfacets_local, sizeof(int));
      for (i = 0; i < nfacets_local; i++) indices[i] = fmap[i].index;

          MPI_Send(indices, nfacets_local, MPI_INT, 0, 98928, mpi->world);
          MPI_Send(local_var, nfacets_local*nvars, MPI_DOUBLE, 0, 98929, mpi->world);

      free(indices);
    }
#endif

    }

  return 0;
}

IBComputeNormalGradient()

int IBComputeNormalGradient	(	void *	s,
		void *	m,
		const double *const	var,
		int	nvars,
		double **const	grad_var
	)

Calculate the normal gradient of a provided variables on the immersed body surface: for each "local facet", i.e., facets (of the immersed body surface) that lie within the local computational domain of this MPI rank, compute the normal gradient.

The variable should be a grid variable with size/layout the same as the solution variable (HyPar::u) with the appropriate number of ghost points.

If the incoming variable has multiple components, this function will calculate normal gradient of each component.

The grad var array should be NULL at input; at output, it will point to an array of size (ImmersedBoundary::nfacets_local X nvars) that contains the normal gradient at each local facet. If there are no local facets, it will remain NULL.

The gradient calculation is currently first-order.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
var	Variable to compute the gradient of
nvars	Number of components in var
grad_var	Array to store the gradient; must be NULL at input

Definition at line 33 of file IBComputeNormalGradient.c.

{
  HyPar             *solver  = (HyPar*)          s;
  MPIVariables      *mpi     = (MPIVariables*)   m; 
  ImmersedBoundary  *IB      = (ImmersedBoundary*) solver->ib;

  if (!solver->flag_ib) return(0);

  if ((*grad_var) != NULL) {
    fprintf(stderr,"Error in IBComputeNormalGradient()\n");
    fprintf(stderr," grad_var is not NULL on rank %d\n",
            mpi->rank );
    return 1;
  }

  int nfacets_local = IB->nfacets_local;
  FacetMap *fmap = IB->fmap;

  if (nfacets_local > 0) {
    (*grad_var) = (double*) calloc (nvars*nfacets_local, sizeof(double));

    for (int n = 0; n < nfacets_local; n++) {
  
      double *alpha;
      int    *nodes, j, k;
  
      double v_c[nvars];
      alpha = &(fmap[n].interp_coeffs[0]);
      nodes = &(fmap[n].interp_nodes[0]);
      _ArraySetValue_(v_c,nvars,0.0);
      for (j=0; j<_IB_NNODES_; j++) {
        for (k=0; k<nvars; k++) {
          v_c[k] += ( alpha[j] * var[nvars*nodes[j]+k] );
        }
      }
  
      double v_ns[nvars];
      alpha = &(fmap[n].interp_coeffs_ns[0]);
      nodes = &(fmap[n].interp_nodes_ns[0]);
      _ArraySetValue_(v_ns,nvars,0.0);
      for (j=0; j<_IB_NNODES_; j++) {
        for (k=0; k<nvars; k++) {
          v_ns[k] += ( alpha[j] * var[nvars*nodes[j]+k] );
        }
      }
  
     double nx = fmap[n].facet->nx;
     double ny = fmap[n].facet->ny;
     double nz = fmap[n].facet->nz;
  
      for (k=0; k<nvars; k++) {
        double dv_dx = (v_ns[k] - v_c[k]) / fmap[n].dx;
        double dv_dy = (v_ns[k] - v_c[k]) / fmap[n].dy;
        double dv_dz = (v_ns[k] - v_c[k]) / fmap[n].dz;
        (*grad_var)[n*nvars+k] = dv_dx*nx + dv_dy*ny + dv_dz*nz;
      }
  
    }

  }

  return(0);
}

IBComputeFacetVar()

int IBComputeFacetVar	(	void *	s,
		void *	m,
		const double *const	var,
		int	nvars,
		double **const	face_var
	)

Calculate the interpolated value of a provided variables on the immersed body surface: for each "local facet", i.e., facets (of the immersed body surface) that lie within the local computational domain of this MPI rank, compute the interpolated value.

The variable should be a grid variable with size/layout the same as the solution variable (HyPar::u) with the appropriate number of ghost points.

If the incoming variable has multiple components, this function will calculate interpolated value of each component.

The grad var array should be NULL at input; at output, it will point to an array of size (ImmersedBoundary::nfacets_local X nvars) that contains the interpolated value at each local facet. If there are no local facets, it will remain NULL.

The interpolation is bi/tri-linear (second-order).

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
var	Variable to compute the interpolated value of
nvars	Number of components in var
face_var	Array to store the interpolated value; must be NULL at input

Definition at line 33 of file IBComputeFacetVar.c.

{
  HyPar             *solver  = (HyPar*)          s;
  MPIVariables      *mpi     = (MPIVariables*)   m; 
  ImmersedBoundary  *IB      = (ImmersedBoundary*) solver->ib;

  if (!solver->flag_ib) return(0);

  if ((*face_var) != NULL) {
    fprintf(stderr,"Error in IBComputeFacetVar()\n");
    fprintf(stderr," face_var is not NULL on rank %d\n",
            mpi->rank );
    return 1;
  }

  int nfacets_local = IB->nfacets_local;
  FacetMap *fmap = IB->fmap;

  if (nfacets_local > 0) {
    (*face_var) = (double*) calloc (nvars*nfacets_local, sizeof(double));

    for (int n = 0; n < nfacets_local; n++) {
  
      double *alpha;
      int    *nodes, j, k;
  
      double *v_c = (*face_var) + n*nvars;
      alpha = &(fmap[n].interp_coeffs[0]);
      nodes = &(fmap[n].interp_nodes[0]);
      _ArraySetValue_(v_c,nvars,0.0);
      for (j=0; j<_IB_NNODES_; j++) {
        for (k=0; k<nvars; k++) {
          v_c[k] += ( alpha[j] * var[nvars*nodes[j]+k] );
        }
      }
    }

  }

  return(0);
}