HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Data Structures | Macros | Functions
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

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

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

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

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

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

#define _IB_NDIMS_   3

Immersed boundaries are implemented only for 3-dimensional simulations.

Definition at line 10 of file immersedboundaries.h.

#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.

#define _IB_XY_   "2d (xy)"

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

Definition at line 16 of file immersedboundaries.h.

#define _IB_XZ_   "2d (xz)"

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

Definition at line 18 of file immersedboundaries.h.

#define _IB_YZ_   "2d (yz)"

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

Definition at line 20 of file immersedboundaries.h.

#define _IB_3D_   "3d"

3D simulation

Definition at line 22 of file immersedboundaries.h.

Function Documentation

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
body3D body to be allocated and read from file
filenameFilename
mMPI object of type MPIVariables
statStatus (0: success; non-0: failure)

Definition at line 21 of file IBReadBodySTL.c.

27 {
28  MPIVariables *mpi = (MPIVariables*) m;
29  int n, ierr;
30  *stat = 0;
31 
32  if ((*body) != NULL) {
33  if (!mpi->rank) {
34  fprintf(stderr,"Error in IBReadBodySTL(): pointer to immersed body not NULL.\n");
35  }
36  *stat = -1;
37  return(0);
38  }
39 
40  /* Rank 0 reads in file */
41  if (!mpi->rank) {
42 
43  FILE *in;
44  in = fopen(filename,"r");
45  if (!in) {
46  fprintf(stderr,"Error in IBReadBodySTL(): file %s not found or cannot be opened.\n",
47  filename);
48  *stat = 1;
49  }
50  else {
51 
52  int nfacets = 0;
53  char word[_MAX_STRING_SIZE_];
54  /* Count number of facets in STL file */
55  ierr = fscanf(in,"%s",word);
56  if (!strcmp(word, "solid")) {
57  nfacets = 0;
58  while (strcmp(word, "endsolid")) {
59  ierr = fscanf(in,"%s",word);
60  if (!strcmp(word, "facet")) nfacets++;
61  }
62  }
63  fclose(in);
64 
65  if (nfacets > 0) {
66  (*body) = (Body3D*) calloc (1,sizeof(Body3D));
67  (*body)->nfacets = nfacets;
68  (*body)->surface = (Facet3D*) calloc (nfacets,sizeof(Facet3D));
69 
70  /* Read in STL data */
71  in = fopen(filename,"r");
72  ierr = fscanf(in,"%s",word);
73  n = 0;
74  while (strcmp(word, "endsolid")) {
75  double t1, t2, t3;
76  ierr = fscanf(in,"%s",word);
77  if (!strcmp(word, "facet")) {
78  ierr = fscanf(in,"%s",word);
79  if (strcmp(word, "normal")) {
80  fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
81  } else {
82  ierr = fscanf(in,"%lf",&t1);
83  ierr = fscanf(in,"%lf",&t2);
84  ierr = fscanf(in,"%lf",&t3);
85  (*body)->surface[n].nx = t1;
86  (*body)->surface[n].ny = t2;
87  (*body)->surface[n].nz = t3;
88  }
89  ierr = fscanf(in,"%s",word);
90  if (strcmp(word, "outer")) {
91  fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
92  *stat = 1;
93  } else {
94  ierr = fscanf(in,"%s",word);
95  if (strcmp(word, "loop")) {
96  fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
97  *stat = 1;
98  } else {
99  ierr = fscanf(in,"%s",word);
100  if (strcmp(word, "vertex")) {
101  fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
102  *stat = 1;
103  } else {
104  ierr = fscanf(in,"%lf",&t1);
105  ierr = fscanf(in,"%lf",&t2);
106  ierr = fscanf(in,"%lf",&t3);
107  (*body)->surface[n].x1 = t1;
108  (*body)->surface[n].y1 = t2;
109  (*body)->surface[n].z1 = t3;
110  }
111  ierr = fscanf(in,"%s",word);
112  if (strcmp(word, "vertex")) {
113  fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
114  *stat = 1;
115  } else {
116  ierr = fscanf(in,"%lf",&t1);
117  ierr = fscanf(in,"%lf",&t2);
118  ierr = fscanf(in,"%lf",&t3);
119  (*body)->surface[n].x2 = t1;
120  (*body)->surface[n].y2 = t2;
121  (*body)->surface[n].z2 = t3;
122  }
123  ierr = fscanf(in,"%s",word);
124  if (strcmp(word, "vertex")) {
125  fprintf(stderr,"Error in IBReadBodySTL(): Illegal keyword read in %s.\n",filename);
126  *stat = 1;
127  } else {
128  ierr = fscanf(in,"%lf",&t1);
129  ierr = fscanf(in,"%lf",&t2);
130  ierr = fscanf(in,"%lf",&t3);
131  (*body)->surface[n].x3 = t1;
132  (*body)->surface[n].y3 = t2;
133  (*body)->surface[n].z3 = t3;
134  }
135  }
136  }
137  n++;
138  }
139  }
140  fclose(in);
141  /* done reading STL file */
142  if (n != nfacets) {
143  fprintf(stderr,"Error in IBReadBodySTL(): Inconsistency in number of facets read.\n");
144  free((*body)->surface);
145  free(*body);
146  *stat = 1;
147  }
148  } else {
149  fprintf(stderr,"Error in IBReadBodySTL(): nfacets = 0!!\n");
150  *stat = 1;
151  }
152 
153  }
154  }
155  MPIBroadcast_integer(stat,1,0,&mpi->world);
156 
157  if ((*stat)) return(0);
158 
159  /* Distribute the body to all MPI ranks */
160  int nfacets;
161  if (!mpi->rank) nfacets = (*body)->nfacets;
162  MPIBroadcast_integer(&nfacets,1,0,&mpi->world);
163 
164  if (mpi->rank) {
165  (*body) = (Body3D*) calloc (1,sizeof(Body3D));
166  (*body)->nfacets = nfacets;
167  (*body)->surface = (Facet3D*) calloc (nfacets,sizeof(Facet3D));
168  }
169 
170  int bufdim = 12;
171  double *buffer = (double*) calloc (bufdim*nfacets,sizeof(double));
172  if (!mpi->rank) {
173  for (n=0; n<nfacets; n++) {
174  buffer[n*bufdim+ 0] = (*body)->surface[n].x1;
175  buffer[n*bufdim+ 1] = (*body)->surface[n].x2;
176  buffer[n*bufdim+ 2] = (*body)->surface[n].x3;
177  buffer[n*bufdim+ 3] = (*body)->surface[n].y1;
178  buffer[n*bufdim+ 4] = (*body)->surface[n].y2;
179  buffer[n*bufdim+ 5] = (*body)->surface[n].y3;
180  buffer[n*bufdim+ 6] = (*body)->surface[n].z1;
181  buffer[n*bufdim+ 7] = (*body)->surface[n].z2;
182  buffer[n*bufdim+ 8] = (*body)->surface[n].z3;
183  buffer[n*bufdim+ 9] = (*body)->surface[n].nx;
184  buffer[n*bufdim+10] = (*body)->surface[n].ny;
185  buffer[n*bufdim+11] = (*body)->surface[n].nz;
186  }
187  }
188  MPIBroadcast_double(buffer,(nfacets*bufdim),0,&mpi->world);
189  if (mpi->rank) {
190  for (n=0; n<nfacets; n++) {
191  (*body)->surface[n].x1 = buffer[n*bufdim+ 0];
192  (*body)->surface[n].x2 = buffer[n*bufdim+ 1];
193  (*body)->surface[n].x3 = buffer[n*bufdim+ 2];
194  (*body)->surface[n].y1 = buffer[n*bufdim+ 3];
195  (*body)->surface[n].y2 = buffer[n*bufdim+ 4];
196  (*body)->surface[n].y3 = buffer[n*bufdim+ 5];
197  (*body)->surface[n].z1 = buffer[n*bufdim+ 6];
198  (*body)->surface[n].z2 = buffer[n*bufdim+ 7];
199  (*body)->surface[n].z3 = buffer[n*bufdim+ 8];
200  (*body)->surface[n].nx = buffer[n*bufdim+ 9];
201  (*body)->surface[n].ny = buffer[n*bufdim+10];
202  (*body)->surface[n].nz = buffer[n*bufdim+11];
203  }
204  }
205  free(buffer);
206 
207  return(0);
208 }
Body3D
Structure defining a body.
Definition: immersedboundaries.h:100
Facet3D
Structure defining a facet.
Definition: immersedboundaries.h:39
MPIBroadcast_double
int MPIBroadcast_double(double *, int, int, void *)
Definition: MPIBroadcast.c:9
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
MPIBroadcast_integer
int MPIBroadcast_integer(int *, int, int, void *)
Definition: MPIBroadcast.c:23
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
body3D body to write to file
filenameFilename
mMPI object of type MPIVariables
rankMPI rank that does the writing
statStatus (0: success; non-0: failure)

Definition at line 23 of file IBWriteBodySTL.c.

30 {
31  MPIVariables *mpi = (MPIVariables*) m;
32 
33  if (mpi->rank == rank) {
34  FILE *out;
35  out = fopen(filename,"w");
36  if (!out) {
37  fprintf(stderr,"Error in IBWriteBodySTL(): Unable to open file ");
38  fprintf(stderr,"%s for writing on rank %d.\n",filename,mpi->rank);
39  *stat = 1;
40  } else {
41  fprintf(out,"solid TEST\n");
42  int n;
43  for (n = 0; n < body->nfacets; n++) {
44  fprintf(out," facet normal %1.16e %1.16e %1.16e\n",
45  body->surface[n].nx,body->surface[n].ny,body->surface[n].nz);
46  fprintf(out," outer loop\n");
47  fprintf(out," vertex %1.16e %1.16e %1.16e\n",
48  body->surface[n].x1,body->surface[n].y1,body->surface[n].z1);
49  fprintf(out," vertex %1.16e %1.16e %1.16e\n",
50  body->surface[n].x2,body->surface[n].y2,body->surface[n].z2);
51  fprintf(out," vertex %1.16e %1.16e %1.16e\n",
52  body->surface[n].x3,body->surface[n].y3,body->surface[n].z3);
53  fprintf(out," endloop\n");
54  fprintf(out," endfacet\n");
55  }
56  fprintf(out,"endsolid TEST\n");
57  *stat = 0;
58  }
59  }
60  return(0);
61 }
Body3D::surface
Facet3D * surface
Definition: immersedboundaries.h:102
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
int IBCleanup ( void *  )

Definition at line 11 of file IBCleanup.c.

12 {
13  ImmersedBoundary *ib = (ImmersedBoundary*) s;
14  if (!ib) return(0);
15 
16  free(ib->body->surface);
17  free(ib->body);
18 
19  if (ib->n_boundary_nodes > 0) free(ib->boundary);
20  if (ib->nfacets_local > 0) free(ib->fmap);
21 
22  return(0);
23 }
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
Body3D::surface
Facet3D * surface
Definition: immersedboundaries.h:102
int IBComputeBoundingBox ( Body3D *  b)

Compute the bounding box for a given body.

Parameters
bThe body

Definition at line 9 of file IBComputeBoundingBox.c.

10 {
11  b->xmin = b->xmax = b->surface[0].x1;
12  b->ymin = b->ymax = b->surface[0].y1;
13  b->zmin = b->zmax = b->surface[0].z1;
14 
15  int n;
16  for (n = 0; n < b->nfacets; n++) {
17  if (b->surface[n].x1 < b->xmin) b->xmin = b->surface[n].x1;
18  if (b->surface[n].x2 < b->xmin) b->xmin = b->surface[n].x2;
19  if (b->surface[n].x3 < b->xmin) b->xmin = b->surface[n].x3;
20 
21  if (b->surface[n].y1 < b->ymin) b->ymin = b->surface[n].y1;
22  if (b->surface[n].y2 < b->ymin) b->ymin = b->surface[n].y2;
23  if (b->surface[n].y3 < b->ymin) b->ymin = b->surface[n].y3;
24 
25  if (b->surface[n].z1 < b->zmin) b->zmin = b->surface[n].z1;
26  if (b->surface[n].z2 < b->zmin) b->zmin = b->surface[n].z2;
27  if (b->surface[n].z3 < b->zmin) b->zmin = b->surface[n].z3;
28 
29  if (b->surface[n].x1 > b->xmax) b->xmax = b->surface[n].x1;
30  if (b->surface[n].x2 > b->xmax) b->xmax = b->surface[n].x2;
31  if (b->surface[n].x3 > b->xmax) b->xmax = b->surface[n].x3;
32 
33  if (b->surface[n].y1 > b->ymax) b->ymax = b->surface[n].y1;
34  if (b->surface[n].y2 > b->ymax) b->ymax = b->surface[n].y2;
35  if (b->surface[n].y3 > b->ymax) b->ymax = b->surface[n].y3;
36 
37  if (b->surface[n].z1 > b->zmax) b->zmax = b->surface[n].z1;
38  if (b->surface[n].z2 > b->zmax) b->zmax = b->surface[n].z2;
39  if (b->surface[n].z3 > b->zmax) b->zmax = b->surface[n].z3;
40  }
41  return(0);
42 }
Body3D::surface
Facet3D * surface
Definition: immersedboundaries.h:102
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
ibImmersed boundary object of type ImmersedBoundary
mMPI object of type MPIVariables
XArray of local spatial coordinates
dimLocal dimensions
ghostsNumber of ghost points

Definition at line 143 of file IBCreateFacetMapping.c.

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::facet
Facet3D * facet
Definition: immersedboundaries.h:68
FacetMap::interp_coeffs
double interp_coeffs[_IB_NNODES_]
Definition: immersedboundaries.h:71
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
Body3D
Structure defining a body.
Definition: immersedboundaries.h:100
Facet3D
Structure defining a facet.
Definition: immersedboundaries.h:39
FacetMap
Structure defining a facet map.
Definition: immersedboundaries.h:67
min3
#define min3(a, b, c)
Definition: math_ops.h:23
isInside
static int isInside(double x, double a, double b)
Definition: IBCreateFacetMapping.c:14
_IB_YZ_
#define _IB_YZ_
Definition: immersedboundaries.h:20
ImmersedBoundary::mode
char mode[_MAX_STRING_SIZE_]
Definition: immersedboundaries.h:164
_IB_3D_
#define _IB_3D_
Definition: immersedboundaries.h:22
Body3D::surface
Facet3D * surface
Definition: immersedboundaries.h:102
_IB_XZ_
#define _IB_XZ_
Definition: immersedboundaries.h:18
FacetMap::interp_coeffs_ns
double interp_coeffs_ns[_IB_NNODES_]
Definition: immersedboundaries.h:73
FacetMap::interp_nodes_ns
int interp_nodes_ns[_IB_NNODES_]
Definition: immersedboundaries.h:72
FacetMap::interp_nodes
int interp_nodes[_IB_NNODES_]
Definition: immersedboundaries.h:70
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
sign
#define sign(a)
Definition: math_ops.h:54
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
_IB_XY_
#define _IB_XY_
Definition: immersedboundaries.h:16
min
#define min(a, b)
Definition: math_ops.h:14
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
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
ibImmersed boundary object of type ImmersedBoundary
mMPI object of type MPIVariables
dim_llocal dimensions
ghostsnumber of ghost points
blankBlanking array: for grid points within the immersed body, this value will be set to 0

Definition at line 158 of file IBIdentifyBoundary.c.

166 {
167  ImmersedBoundary *IB = (ImmersedBoundary*) ib;
168  MPIVariables *mpi = (MPIVariables*) m;
169  Body3D *body = IB->body;
170 
171  int imax = dim_l[0],
172  jmax = dim_l[1],
173  kmax = dim_l[2];
174 
175  int n_boundary_nodes = CountBoundaryPoints(imax,jmax,kmax,ghosts,blank);
176  IB->n_boundary_nodes = n_boundary_nodes;
177  if (n_boundary_nodes == 0) IB->boundary = NULL;
178  else {
179  IB->boundary = (IBNode*) calloc (n_boundary_nodes, sizeof(IBNode));
180  int check = SetBoundaryPoints(imax,jmax,kmax,ghosts,blank,IB->boundary);
181  if (check != n_boundary_nodes) {
182  fprintf(stderr,"Error in IBIdentifyBoundary(): Inconsistency encountered when setting boundary indices. ");
183  fprintf(stderr,"on rank %d.\n",mpi->rank);
184  fprintf(stderr,"SetBoundaryPoints() returned %d, while n_boundary_nodes is %d.\n",check,n_boundary_nodes);
185  }
186  }
187 
188  return(n_boundary_nodes);
189 }
IBNode
Structure defining an immersed boundary node.
Definition: immersedboundaries.h:125
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
Body3D
Structure defining a body.
Definition: immersedboundaries.h:100
CountBoundaryPoints
static int CountBoundaryPoints(int imax, int jmax, int kmax, int ghosts, double *blank)
Definition: IBIdentifyBoundary.c:15
SetBoundaryPoints
static int SetBoundaryPoints(int imax, int jmax, int kmax, int ghosts, double *blank, void *b)
Definition: IBIdentifyBoundary.c:80
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
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
XArray of global spatial coordinates
dimglobal dimensions
ibImmersed boundary object of type ImmersedBoundary

Definition at line 24 of file IBIdentifyMode.c.

29 {
30  ImmersedBoundary *IB = (ImmersedBoundary*) ib;
31  Body3D *body = IB->body;
32 
33  double *x = X, *y = x + dim[0], *z = y + dim[1];
34 
35  double grid_xmin = x[0],
36  grid_xmax = x[dim[0]-1],
37  grid_ymin = y[0],
38  grid_ymax = y[dim[1]-1],
39  grid_zmin = z[0],
40  grid_zmax = z[dim[2]-1];
41 
42  if ( (grid_xmin > body->xmin) && (grid_xmax < body->xmax) ) strcpy(IB->mode,_IB_YZ_);
43  else if ( (grid_ymin > body->ymin) && (grid_ymax < body->ymax) ) strcpy(IB->mode,_IB_XZ_);
44  else if ( (grid_zmin > body->zmin) && (grid_zmax < body->zmax) ) strcpy(IB->mode,_IB_XY_);
45  else strcpy(IB->mode,_IB_3D_);
46 
47  return(0);
48 }
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
Body3D
Structure defining a body.
Definition: immersedboundaries.h:100
_IB_YZ_
#define _IB_YZ_
Definition: immersedboundaries.h:20
ImmersedBoundary::mode
char mode[_MAX_STRING_SIZE_]
Definition: immersedboundaries.h:164
_IB_3D_
#define _IB_3D_
Definition: immersedboundaries.h:22
_IB_XZ_
#define _IB_XZ_
Definition: immersedboundaries.h:18
_IB_XY_
#define _IB_XY_
Definition: immersedboundaries.h:16
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
ibImmersed boundary object of type ImmersedBoundary
mMPI object of type MPIVariables
XArray of (local) spatial coordinates
large_distanceA large distance
dim_lInteger array of local grid size in each spatial dimension
ghostsNumber of ghost points

Definition at line 25 of file IBNearestFacetNormal.c.

33 {
34  ImmersedBoundary *IB = (ImmersedBoundary*) ib;
35  MPIVariables *mpi = (MPIVariables*) m;
36  Body3D *body = IB->body;
37  Facet3D *surface = body->surface;
38  IBNode *boundary = IB->boundary;
39 
40  double eps = IB->tolerance;
41  int nb = IB->n_boundary_nodes,
42  nf = body->nfacets;
43 
44  int i, j, k, dg, n;
45  for (dg = 0; dg < nb; dg++) {
46  i = boundary[dg].i;
47  j = boundary[dg].j;
48  k = boundary[dg].k;
49 
50  double xp, yp, zp;
51  _GetCoordinate_(0,i,dim_l,ghosts,X,xp);
52  _GetCoordinate_(1,j,dim_l,ghosts,X,yp);
53  _GetCoordinate_(2,k,dim_l,ghosts,X,zp);
54  boundary[dg].x = xp;
55  boundary[dg].y = yp;
56  boundary[dg].z = zp;
57 
58  double dist_min = large_distance;
59  int n_min = -1;
60 
61  for (n = 0; n < nf; n++) {
62  double x1, x2, x3;
63  double y1, y2, y3;
64  double z1, z2, z3;
65  x1 = surface[n].x1; x2 = surface[n].x2; x3 = surface[n].x3;
66  y1 = surface[n].y1; y2 = surface[n].y2; y3 = surface[n].y3;
67  z1 = surface[n].z1; z2 = surface[n].z2; z3 = surface[n].z3;
68  double dist = surface[n].nx*(xp-surface[n].x1)
69  + surface[n].ny*(yp-surface[n].y1)
70  + surface[n].nz*(zp-surface[n].z1);
71  if (dist > 0) continue;
72  if (absolute(dist) < dist_min) {
73  short is_it_in = 0;
74  double x_int, y_int, z_int;
75  x_int = xp - dist * surface[n].nx;
76  y_int = yp - dist * surface[n].ny;
77  z_int = zp - dist * surface[n].nz;
78  if (absolute(surface[n].nx) > eps) {
79  double den = (z2-z3)*(y1-y3)-(y2-y3)*(z1-z3);
80  double l1, l2, l3;
81  l1 = ((y2-y3)*(z3-z_int)-(z2-z3)*(y3-y_int)) / den;
82  l2 = ((z1-z3)*(y3-y_int)-(y1-y3)*(z3-z_int)) / den;
83  l3 = 1 - l1 - l2;
84  if ((l1 > -eps) && (l2 > -eps) && (l3 > -eps)) is_it_in = 1;
85  } else if (absolute(surface[n].ny) > eps) {
86  double den = (x2-x3)*(z1-z3)-(z2-z3)*(x1-x3);
87  double l1, l2, l3;
88  l1 = ((z2-z3)*(x3-x_int)-(x2-x3)*(z3-z_int)) / den;
89  l2 = ((x1-x3)*(z3-z_int)-(z1-z3)*(x3-x_int)) / den;
90  l3 = 1 - l1 - l2;
91  if ((l1 > -eps) && (l2 > -eps) && (l3 > -eps)) is_it_in = 1;
92  } else {
93  double den = (y2-y3)*(x1-x3)-(x2-x3)*(y1-y3);
94  double l1, l2, l3;
95  l1 = ((x2-x3)*(y3-y_int)-(y2-y3)*(x3-x_int)) / den;
96  l2 = ((y1-y3)*(x3-x_int)-(x1-x3)*(y3-y_int)) / den;
97  l3 = 1 - l1 - l2;
98  if ((l1 > -eps) && (l2 > -eps) && (l3 > -eps)) is_it_in = 1;
99  }
100  if (is_it_in) {
101  dist_min = absolute(dist);
102  n_min = n;
103  }
104  }
105  }
106  if (n_min == -1) {
107  for (n = 0; n < nf; n++) {
108  double dist = surface[n].nx*(xp-surface[n].x1)
109  + surface[n].ny*(yp-surface[n].y1)
110  + surface[n].nz*(zp-surface[n].z1);
111  if (dist > eps) continue;
112  else {
113  if (absolute(dist) < dist_min) {
114  dist_min = absolute(dist);
115  n_min = n;
116  }
117  }
118  }
119  }
120 
121  if (n_min == -1) {
122  fprintf(stderr,"Error in IBNearestFacetNormal(): no nearest normal found for boundary node (%d,%d,%d) ",i,j,k);
123  fprintf(stderr,"on rank %d.\n",mpi->rank);
124  return(1);
125  } else boundary[dg].face = &surface[n_min];
126  }
127 
128  return(0);
129 }
IBNode
Structure defining an immersed boundary node.
Definition: immersedboundaries.h:125
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
Body3D
Structure defining a body.
Definition: immersedboundaries.h:100
Facet3D
Structure defining a facet.
Definition: immersedboundaries.h:39
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
Body3D::surface
Facet3D * surface
Definition: immersedboundaries.h:102
absolute
#define absolute(a)
Definition: math_ops.h:32
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
IBNode::face
Facet3D * face
Definition: immersedboundaries.h:133
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
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
mMPI object of type MPIVariables
ibImmersed boundary object of type ImmersedBoundary
local_varLocal array
global_varArray to store the gradient; must be NULL at input
nvarsNumber of components in var

Definition at line 22 of file IBAssembleGlobalFacetData.c.

28 {
29  MPIVariables *mpi = (MPIVariables*) m;
30  ImmersedBoundary *IB = (ImmersedBoundary*) ib;
31 
32  if ((*global_var) != NULL) {
33  fprintf(stderr,"Error in IBAssembleGlobalFacetData()\n");
34  fprintf(stderr," global_var is not NULL on rank %d\n",
35  mpi->rank );
36  return 1;
37  }
38 
39  int nfacets_local = IB->nfacets_local;
40  FacetMap *fmap = IB->fmap;
41 
42  if ((nfacets_local == 0) && (local_var != NULL)) {
43  fprintf(stderr, "Error in IBAssembleGlobalFacetData()\n");
44  fprintf(stderr, " nfacets_local is 0 but local_var is not NULL!\n");
45  return 1;
46  }
47  if ((nfacets_local > 0) && (local_var == NULL)) {
48  fprintf(stderr, "Error in IBAssembleGlobalFacetData()\n");
49  fprintf(stderr, " nfacets_local > 0 but local_var is NULL!\n");
50  return 1;
51  }
52 
53 #ifndef serial
54  MPI_Status status;
55 #endif
56 
57  if (!mpi->rank) {
58 
59  int nfacets_global = IB->body->nfacets;
60 
61  /* allocate arrays for whole domain */
62  *global_var = (double*) calloc (nfacets_global*nvars, sizeof(double));
63  _ArraySetValue_((*global_var), (nfacets_global*nvars), 0.0);
64 
65  /* check array - to avoid double counting of facets */
66  int *check = (int*) calloc (nfacets_global, sizeof(int));
67  _ArraySetValue_(check, nfacets_global, 0);
68 
69  /* local data */
70  for (int n = 0; n < nfacets_local; n++) {
71  _ArrayAXPY_((local_var+n), 1.0, ((*global_var)+fmap[n].index), nvars);
72  check[fmap[n].index]++;
73  }
74 
75 #ifndef serial
76  for (int proc = 1; proc < mpi->nproc; proc++) {
77 
78  int nf_incoming;
79  MPI_Recv(&nf_incoming, 1, MPI_INT, proc, 98927, MPI_COMM_WORLD, &status);
80 
81  if (nf_incoming > 0) {
82 
83  int *indices_incoming = (int*) calloc(nf_incoming, sizeof(int));
84  double *var_incoming = (double*) calloc(nf_incoming*nvars, sizeof(double));
85 
86  MPI_Recv(indices_incoming, nf_incoming, MPI_INT, proc, 98928, mpi->world, &status);
87  MPI_Recv(var_incoming, nf_incoming*nvars, MPI_DOUBLE, proc, 98929, mpi->world, &status);
88 
89  for (int n = 0; n < nf_incoming; n++) {
90  _ArrayAXPY_((var_incoming+n), 1.0, ((*global_var)+indices_incoming[n]), nvars);
91  check[indices_incoming[n]]++;
92  }
93 
94  free(indices_incoming);
95  free(var_incoming);
96  }
97  }
98 #endif
99 
100  for (int n = 0; n < nfacets_global; n++) {
101  if (check[n] == 0) {
102  fprintf(stderr, "Error in IBAssembleGlobalFacetData()\n");
103  fprintf(stderr, " No data received for facet %d\n", n);
104  return 1;
105  } else {
106  _ArrayScale1D_(((*global_var)+n), (1.0/(double)check[n]), nvars);
107  }
108  }
109 
110  } else {
111 
112 #ifndef serial
113 
114  MPI_Send(&nfacets_local, 1, MPI_INT, 0, 98927, MPI_COMM_WORLD);
115 
116  if (nfacets_local > 0) {
117 
118  int i, *indices = (int*) calloc (nfacets_local, sizeof(int));
119  for (i = 0; i < nfacets_local; i++) indices[i] = fmap[i].index;
120 
121  MPI_Send(indices, nfacets_local, MPI_INT, 0, 98928, mpi->world);
122  MPI_Send(local_var, nfacets_local*nvars, MPI_DOUBLE, 0, 98929, mpi->world);
123 
124  free(indices);
125  }
126 #endif
127 
128  }
129 
130  return 0;
131 }
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
FacetMap
Structure defining a facet map.
Definition: immersedboundaries.h:67
_ArrayScale1D_
#define _ArrayScale1D_(x, a, size)
Definition: arrayfunctions.h:190
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
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
sSolver object of type HyPar
mMPI object of type MPIVariables
varVariable to compute the gradient of
nvarsNumber of components in var
grad_varArray to store the gradient; must be NULL at input

Definition at line 33 of file IBComputeNormalGradient.c.

39 {
40  HyPar *solver = (HyPar*) s;
41  MPIVariables *mpi = (MPIVariables*) m;
42  ImmersedBoundary *IB = (ImmersedBoundary*) solver->ib;
43 
44  if (!solver->flag_ib) return(0);
45 
46  if ((*grad_var) != NULL) {
47  fprintf(stderr,"Error in IBComputeNormalGradient()\n");
48  fprintf(stderr," grad_var is not NULL on rank %d\n",
49  mpi->rank );
50  return 1;
51  }
52 
53  int nfacets_local = IB->nfacets_local;
54  FacetMap *fmap = IB->fmap;
55 
56  if (nfacets_local > 0) {
57  (*grad_var) = (double*) calloc (nvars*nfacets_local, sizeof(double));
58 
59  for (int n = 0; n < nfacets_local; n++) {
60 
61  double *alpha;
62  int *nodes, j, k;
63 
64  double v_c[nvars];
65  alpha = &(fmap[n].interp_coeffs[0]);
66  nodes = &(fmap[n].interp_nodes[0]);
67  _ArraySetValue_(v_c,nvars,0.0);
68  for (j=0; j<_IB_NNODES_; j++) {
69  for (k=0; k<nvars; k++) {
70  v_c[k] += ( alpha[j] * var[nvars*nodes[j]+k] );
71  }
72  }
73 
74  double v_ns[nvars];
75  alpha = &(fmap[n].interp_coeffs_ns[0]);
76  nodes = &(fmap[n].interp_nodes_ns[0]);
77  _ArraySetValue_(v_ns,nvars,0.0);
78  for (j=0; j<_IB_NNODES_; j++) {
79  for (k=0; k<nvars; k++) {
80  v_ns[k] += ( alpha[j] * var[nvars*nodes[j]+k] );
81  }
82  }
83 
84  double nx = fmap[n].facet->nx;
85  double ny = fmap[n].facet->ny;
86  double nz = fmap[n].facet->nz;
87 
88  for (k=0; k<nvars; k++) {
89  double dv_dx = (v_ns[k] - v_c[k]) / fmap[n].dx;
90  double dv_dy = (v_ns[k] - v_c[k]) / fmap[n].dy;
91  double dv_dz = (v_ns[k] - v_c[k]) / fmap[n].dz;
92  (*grad_var)[n*nvars+k] = dv_dx*nx + dv_dy*ny + dv_dz*nz;
93  }
94 
95  }
96 
97  }
98 
99  return(0);
100 }
FacetMap::facet
Facet3D * facet
Definition: immersedboundaries.h:68
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
FacetMap::interp_coeffs
double interp_coeffs[_IB_NNODES_]
Definition: immersedboundaries.h:71
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
FacetMap
Structure defining a facet map.
Definition: immersedboundaries.h:67
HyPar::ib
void * ib
Definition: hypar.h:443
HyPar::flag_ib
int flag_ib
Definition: hypar.h:441
FacetMap::interp_coeffs_ns
double interp_coeffs_ns[_IB_NNODES_]
Definition: immersedboundaries.h:73
FacetMap::interp_nodes_ns
int interp_nodes_ns[_IB_NNODES_]
Definition: immersedboundaries.h:72
FacetMap::interp_nodes
int interp_nodes[_IB_NNODES_]
Definition: immersedboundaries.h:70
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_IB_NNODES_
#define _IB_NNODES_
Definition: immersedboundaries.h:13
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
sSolver object of type HyPar
mMPI object of type MPIVariables
varVariable to compute the interpolated value of
nvarsNumber of components in var
face_varArray to store the interpolated value; must be NULL at input

Definition at line 33 of file IBComputeFacetVar.c.

40 {
41  HyPar *solver = (HyPar*) s;
42  MPIVariables *mpi = (MPIVariables*) m;
43  ImmersedBoundary *IB = (ImmersedBoundary*) solver->ib;
44 
45  if (!solver->flag_ib) return(0);
46 
47  if ((*face_var) != NULL) {
48  fprintf(stderr,"Error in IBComputeFacetVar()\n");
49  fprintf(stderr," face_var is not NULL on rank %d\n",
50  mpi->rank );
51  return 1;
52  }
53 
54  int nfacets_local = IB->nfacets_local;
55  FacetMap *fmap = IB->fmap;
56 
57  if (nfacets_local > 0) {
58  (*face_var) = (double*) calloc (nvars*nfacets_local, sizeof(double));
59 
60  for (int n = 0; n < nfacets_local; n++) {
61 
62  double *alpha;
63  int *nodes, j, k;
64 
65  double *v_c = (*face_var) + n*nvars;
66  alpha = &(fmap[n].interp_coeffs[0]);
67  nodes = &(fmap[n].interp_nodes[0]);
68  _ArraySetValue_(v_c,nvars,0.0);
69  for (j=0; j<_IB_NNODES_; j++) {
70  for (k=0; k<nvars; k++) {
71  v_c[k] += ( alpha[j] * var[nvars*nodes[j]+k] );
72  }
73  }
74  }
75 
76  }
77 
78  return(0);
79 }
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
FacetMap::interp_coeffs
double interp_coeffs[_IB_NNODES_]
Definition: immersedboundaries.h:71
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
FacetMap
Structure defining a facet map.
Definition: immersedboundaries.h:67
HyPar::ib
void * ib
Definition: hypar.h:443
HyPar::flag_ib
int flag_ib
Definition: hypar.h:441
FacetMap::interp_nodes
int interp_nodes[_IB_NNODES_]
Definition: immersedboundaries.h:70
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_IB_NNODES_
#define _IB_NNODES_
Definition: immersedboundaries.h:13