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

Find the nearest facet for immersed boundary points. More...

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

Go to the source code of this file.

Functions

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

Detailed Description

Find the nearest facet for immersed boundary points.

Author
Debojyoti Ghosh

Definition in file IBNearestFacetNormal.c.

Function Documentation

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 {
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 }
Structure defining an immersed boundary node.
Structure containing variables for immersed boundary implementation.
Structure defining a body.
Structure defining a facet.
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
Facet3D * surface
#define absolute(a)
Definition: math_ops.h:32
Structure of MPI-related variables.
Facet3D * face