embed.h

Go to the documentation of this file.

/** @file embed.h
 */
/**
# Embedded boundaries
 
Boundaries of general shape can be described using an integral
(i.e. finite volume) formulation which takes into account the volume
and area fractions of intersection of the embedded boundary with the
Cartesian mesh. 
 
We will need to deal with volume fractions. Interpolations (for
Dirichlet boundary conditions) assume a 5x5 stencil. */
 
#include "fractions.h"
#define BGHOSTS 2
#define EMBED 1
 
/**
The volume and area fractions are stored in these fields. */
 
scalar cs[];
vector fs[];
 
double (* metric_embed_factor) (Point, coord) = NULL;
 
/**
Embedded boundary operators specific to trees are defined in this
file. */
 
#if TREE
# include "embed-tree.h"
#endif
 
/**
## Operator overloading
 
Several standard operators, defined in [common.h]() need to be tuned
to take into account the embedded fractions. 
 
The *SEPS* constant is used to avoid division by zero. */
 
#undef SEPS
#define SEPS 1e-30
 
/**
When combining third-order Dirichlet conditions and [approximate
projections](navier-stokes/centered.h#approximate-projection),
instabilities may occur due to feedback between the pressure mode and
the velocity, amplified by the third-order derivative. This can be
stabilised using the weighted average below when computing face
velocities. The corresponding test case is [test/uf.c](). Note that if
only second-order Dirichlet fluxes are used, simple averaging is
stable. */
 
#define cs_avg(a,i,j,k)             \
  ((a[i,j,k]*(1.5 + cs[i,j,k]) + a[i-1,j,k]*(1.5 + cs[i-1,j,k]))/ \
   (cs[i,j,k] + cs[i-1,j,k] + 3.))
 
/**
Face gradients and face values, computed from cell-centered values
must be tuned to take into account the area fractions of the embedded
boundary. This follows the procedure described in [Johansen and
Colella, 1998](#johansen1998), figure 3 and equation (16) in
particular. Note that this is only used when using second-order
fluxes. */
 
#if dimension == 2
#define face_condition(fs, cs)            \
  (fs.x[i,j] > 0.5 && fs.y[i,j + (j < 0)] && fs.y[i-1,j + (j < 0)] && \
   cs[i,j] && cs[i-1,j])
 
for (int _d = 0; _d < dimension; _d++)
static inline double embed_face_gradient_x (Point point, scalar a, int i)
{
  int j = sign(fs.x[i,1] - fs.x[i,-1]);
  assert (cs[i] && cs[i-1]);
  if (face_condition (fs, cs))
    return ((1. + fs.x[i])*(a[i] - a[i-1]) +
      (1. - fs.x[i])*(a[i,j] - a[i-1,j]))/(2.*Delta);
  return (a[i] - a[i-1])/Delta;
}
 
for (int _d = 0; _d < dimension; _d++)
static inline double embed_face_value_x (Point point, scalar a, int i)
{
  int j = sign(fs.x[i,1] - fs.x[i,-1]);
  return face_condition (fs, cs) ?
    ((1. + fs.x[i])*cs_avg(a,i,0,0) + (1. - fs.x[i])*cs_avg(a,i,j,0))/2. :
    cs_avg(a,i,0,0);
}
 
/**
The generalisation to 3D is a bit more complicated. See Fig. 1 of
[Schwartz et al, 2006](#schwartz2006). */
 
#else // dimension == 3
for (int _d = 0; _d < dimension; _d++)
static inline coord embed_face_barycentre_z (Point point, int i)
{
  // Young's normal calculation
  coord n1 = {0};
  double nn = 0.;
  scalar f = fs.z;
  foreach_dimension(2) {
    n1.x = (f[-1,-1,i] + 2.*f[-1,0,i] + f[-1,1,i] -
      f[+1,-1,i] - 2.*f[+1,0,i] - f[+1,1,i]);
    nn += fabs(n1.x);
  }
  if (!nn)
    return (coord){0.,0.,0.};
  foreach_dimension(2)
    n1.x /= nn;
  // Position `p` of the face barycentre
  coord n, p1, p;
  ((double *)&n)[0] = n1.x, ((double *)&n)[1] = n1.y;
  double alpha = line_alpha (f[0,0,i], n);
  line_center (n, alpha, f[0,0,i], &p1);
  p.x = ((double *)&p1)[0], p.y = ((double *)&p1)[1], p.z = 0.;
  return p;
}
 
/**
The macro below defines the condition which must be verified when
applying the second-order flux interpolation and/or face
averaging. The goal is to avoid using values from areas of the mesh
which are not topologically connected. Doing so would couple
disconnected subproblems (for example when solving a Laplacian) which
would most probably lead to lack of convergence. This is very
important for robustness when dealing with complex boundaries. */
 
#define face_condition(fs, cs)            \
  (fs.x[i,j,k] > 0.5 && (fs.x[i,j,0] > 0.5 || fs.x[i,0,k] > 0.5) && \
   fs.y[i,j + (j < 0),0] && fs.y[i-1,j + (j < 0),0] &&      \
   fs.y[i,j + (j < 0),k] && fs.y[i-1,j + (j < 0),k] &&      \
   fs.z[i,0,k + (k < 0)] && fs.z[i-1,0,k + (k < 0)] &&      \
   fs.z[i,j,k + (k < 0)] && fs.z[i-1,j,k + (k < 0)] &&      \
   cs[i-1,j,0] && cs[i-1,0,k] && cs[i-1,j,k] &&       \
   cs[i,j,0] && cs[i,0,k] && cs[i,j,k])
 
for (int _d = 0; _d < dimension; _d++)
static inline double embed_face_gradient_x (Point point, scalar a, int i)
{
  assert (cs[i] && cs[i-1]);
  coord p = embed_face_barycentre_x (point, i);
  // Bilinear interpolation of the gradient (see Fig. 1 of Schwartz et al., 2006)
  int j = sign(p.y), k = sign(p.z);
  if (face_condition(fs, cs)) {
    p.y = fabs(p.y), p.z = fabs(p.z);
    return (((a[i,0,0] - a[i-1,0,0])*(1. - p.y) +
       (a[i,j,0] - a[i-1,j,0])*p.y)*(1. - p.z) + 
      ((a[i,0,k] - a[i-1,0,k])*(1. - p.y) +
       (a[i,j,k] - a[i-1,j,k])*p.y)*p.z)/Delta;
  }
  return (a[i] - a[i-1])/Delta;
}
 
for (int _d = 0; _d < dimension; _d++)
static inline double embed_face_value_x (Point point, scalar a, int i)
{
  coord p = embed_face_barycentre_x (point, i);
  // Bilinear interpolation
  int j = sign(p.y), k = sign(p.z);
  if (face_condition(fs, cs)) {
    p.y = fabs(p.y), p.z = fabs(p.z);
    return ((cs_avg(a,i,0,0)*(1. - p.y) + cs_avg(a,i,j,0)*p.y)*(1. - p.z) + 
      (cs_avg(a,i,0,k)*(1. - p.y) + cs_avg(a,i,j,k)*p.y)*p.z);
  }
  return cs_avg(a,i,0,0);
}
#endif // dimension == 3
 
/**
We use the functions above to redefine the face gradient macros. Note
that the second-order face gradients and averaging are used only if
the corresponding scalar attribute below (`third` because of
third-order accuracy when using Dirichlet conditions, see
[test/neumann.c]) is set to `true`. The default is `false`. */
 
attribute {
  bool third;
}
 
#undef face_gradient_x
#define face_gradient_x(a,i)          \
  (a.third && fs.x[i] < 1. && fs.x[i] > 0. ?      \
   embed_face_gradient_x (point, a, i) :      \
   (a[i] - a[i-1])/Delta)
 
#undef face_gradient_y
#define face_gradient_y(a,i)          \
  (a.third && fs.y[0,i] < 1. && fs.y[0,i] > 0. ?    \
   embed_face_gradient_y (point, a, i) :      \
   (a[0,i] - a[0,i-1])/Delta)
 
#undef face_gradient_z
#define face_gradient_z(a,i)          \
  (a.third && fs.z[0,0,i] < 1. && fs.z[0,0,i] > 0. ?    \
   embed_face_gradient_z (point, a, i) :      \
   (a[0,0,i] - a[0,0,i-1])/Delta)
 
#undef face_value
#define face_value(a,i)             \
  (a.third && fs.x[i] < 1. && fs.x[i] > 0. ?        \
   embed_face_value_x (point, a, i) :         \
   cs_avg(a,i,0,0))
 
/**
The centered gradient must not use values of fields entirely contained
within the embedded boundary (for which *cs* is zero). */
 
#undef center_gradient
#define center_gradient(a) (fs.x[] && fs.x[1] ? (a[1] - a[-1])/(2.*Delta) : \
          fs.x[1] ? (a[1] - a[])/Delta :        \
          fs.x[]  ? (a[] - a[-1])/Delta : 0.)
 
/**
## Utility functions for the geometry of embedded boundaries
 
For a cell containing a fragment of embedded boundary (i.e. for which
\f$0 < cs < 1\f$), *embed_geometry()* returns the area of the fragment,
the relative position *p* of the barycenter of the fragment and the
boundary normal *n*. */
 
static inline
double embed_geometry (Point point, coord * p, coord * n)
{
  *n = facet_normal (point, cs, fs);
  double alpha = plane_alpha (cs[], *n);
  double area = plane_area_center (*n, alpha, p);
  normalize (n);
  return area;
}
 
/**
This function and the macro below shift the position $(x1,y1,z1)$ to
the position of the barycenter of the embedded fragment.  */
 
static inline
double embed_area_center (Point point, double * x1, double * y1, double * z1)
{
  double area = 0.;
  if (cs[] > 0. && cs[] < 1.) {
    coord n, p;
    area = embed_geometry (point, &p, &n);
    *x1 += p.x*Delta, *y1 += p.y*Delta, *z1 += p.z*Delta;
  }
  return area;
}
 
#define embed_pos() embed_area_center (point, &x, &y, &z)
 
/**
This function returns the value of field *s* interpolated linearly at
the barycenter *p* of the fragment of embedded boundary contained
within the cell. */
 
double embed_interpolate (Point point, scalar s, coord p)
{
  assert (dimension == 2);
  int i = sign(p.x), j = sign(p.y);
  if (cs[i] && cs[0,j] && cs[i,j])
    // bilinear interpolation when all neighbors are defined
    return ((s[]*(1. - fabs(p.x)) + s[i]*fabs(p.x))*(1. - fabs(p.y)) + 
      (s[0,j]*(1. - fabs(p.x)) + s[i,j]*fabs(p.x))*fabs(p.y));
  else {
    // linear interpolation with gradients biased toward the
    // cells which are defined
    double val = s[];
    for (int _d = 0; _d < dimension; _d++) {
      int i = sign(p.x);
      if (cs[i])
  val += fabs(p.x)*(s[i] - s[]);
      else if (cs[-i])
  val += fabs(p.x)*(s[] - s[-i]);
    }
    return val;
  }
}
 
/**
This function "removes" (by setting their volume fraction to zero)
cells which have inconsistent volume/surface fractions. This is
important to guarantee the robustness of the solution for complex (and
under-resolved) boundaries. The functions returns the number of cells
removed. */
 
struct Cleanup {
  scalar c;
  vector s;
  double smin;   // minimum surface fraction (optional)
  bool opposite; // whether to eliminate 'thin tubes' (optional)
};
 
trace
int fractions_cleanup (scalar c, vector s,
           double smin = 0., bool opposite = false)
{
  
  /**
  Since both surface and volume fractions are altered, iterations are
  needed. This reflects the fact that changes are coupled spatially
  through the topology of the domain: for examples, long, unresolved
  "filaments" may need many iterations to be fully removed. */
  
  int changed = 1, schanged = 0, i;
  for (i = 0; i < 100 && changed; i++) {
 
    /**
    Face fractions of empty cells must be zero. */
   
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      if (s.x[] && ((!c[] || !c[-1]) || s.x[] < smin))
  s.x[] = 0.;
 
    changed = 0;
    foreach(reduction(+:changed))
      if (c[] > 0. && c[] < 1.) {
  int n = 0;
  for (int _d = 0; _d < dimension; _d++) {
    for (int i = 0; i <= 1; i++)
      if (s.x[i] > 0.)
        n++;
 
    /**
    If opposite surface fractions are zero (and volume fraction
    is non-zero), then we are dealing with a thin "tube", which
    we just remove because it can sometimes lead to
    non-convergence when
    [projecting](navier-stokes/centered.h#approximate-projection)
    the velocity field. */
 
    if (opposite && s.x[] == 0. && s.x[1] == 0.)
      c[] = 0., changed++;
  }
 
  /**
  The number of "non-empty" faces (i.e. faces which have a
  surface fraction larger than epsilon) cannot be smaller than
  the dimension (the limiting cases correspond to a triangle in
  2D and a tetrahedron in 3D). */
  
  if (n < dimension)
    c[] = 0., changed++;
      }
 
    schanged += changed;
  }
  if (changed)
    fprintf (stderr, "src/embed.h:%d: warning: fractions_cleanup() did not converge after "
       "%d iterations\n", LINENO, i);
  return schanged;
}
  
/**
## Dirichlet boundary condition
 
This function returns the gradient of scalar *s*, normal to the
embedded boundary defined by *cs*, of unit normal vector *n*
(normalised using the Euclidean norm, not the box norm) and of
centroid *p*. The Dirichlet boundary condition *bc* is imposed on the
embedded boundary.
 
The calculation follows [Johansen and Colella, 1998](#johansen1998)
and is summarised in the figure below (see also Figure 4 of Johansen
and Colella and Figure 2 of [Schwartz et al, 2006](#schwartz2006) for
the 3D implementation).
 
![Third-order normal gradient scheme](figures/dirichlet_gradient.svg) 
 
For degenerate cases, a non-zero value of *coef* is returned and
`coef*s[]` must be added to the value returned to obtain the gradient. */
 
#define quadratic(x,a1,a2,a3) \
  (((a1)*((x) - 1.) + (a3)*((x) + 1.))*(x)/2. - (a2)*((x) - 1.)*((x) + 1.))
 
for (int _d = 0; _d < dimension; _d++)
static inline double dirichlet_gradient_x (Point point, scalar s, scalar cs,
             coord n, coord p, double bc,
             double * coef)
{
  for (int _d = 0; _d < dimension; _d++)
    n.x = - n.x;
  double d[2], v[2] = {nodata,nodata};
  bool defined = true;
  for (int _d = 0; _d < dimension; _d++)
    if (defined && !fs.x[(n.x > 0.)])
      defined = false;
  if (defined)
    for (int l = 0; l <= 1; l++) {
      int i = (l + 1)*sign(n.x);
      d[l] = (i - p.x)/n.x;
      double y1 = p.y + d[l]*n.y;
      int j = y1 > 0.5 ? 1 : y1 < -0.5 ? -1 : 0;
      y1 -= j;
#if dimension == 2
      if (fs.x[i + (i < 0),j] && fs.y[i,j] && fs.y[i,j+1] &&
    cs[i,j-1] && cs[i,j] && cs[i,j+1])
  v[l] = quadratic (y1, (s[i,j-1]), (s[i,j]), (s[i,j+1]));
#else // dimension == 3
      double z = p.z + d[l]*n.z;
      int k = z > 0.5 ? 1 : z < -0.5 ? -1 : 0;
      z -= k;
      bool defined = fs.x[i + (i < 0),j,k];
      for (int m = -1; m <= 1 && defined; m++)
  if (!fs.y[i,j,k+m] || !fs.y[i,j+1,k+m] ||
      !fs.z[i,j+m,k] || !fs.z[i,j+m,k+1] ||
      !cs[i,j+m,k-1] || !cs[i,j+m,k] || !cs[i,j+m,k+1])
    defined = false;
      if (defined)
  // bi-quadratic interpolation
  v[l] =
    quadratic (z,
         quadratic (y1,
        (s[i,j-1,k-1]), (s[i,j,k-1]), (s[i,j+1,k-1])),
         quadratic (y1,
        (s[i,j-1,k]),   (s[i,j,k]),   (s[i,j+1,k])),
         quadratic (y1,
        (s[i,j-1,k+1]), (s[i,j,k+1]), (s[i,j+1,k+1])));
#endif // dimension == 3
      else
  break;
    }
  if (v[0] == nodata) {
 
    /**
    This is a degenerate case, we use the boundary value and the
    cell-center value to define the gradient. */
  
    d[0] = max(1e-3, fabs(p.x/n.x));
    *coef = - 1./(d[0]*Delta);
    return bc/(d[0]*Delta);
  }
 
  /**
  For non-degenerate cases, the gradient is obtained using either
  second- or third-order estimates. */
  
  *coef = 0.;
  if (v[1] != nodata) // third-order gradient
    return (d[1]*(bc - v[0])/d[0] - d[0]*(bc - v[1])/d[1])/((d[1] - d[0])*Delta);
  return (bc - v[0])/(d[0]*Delta); // second-order gradient
}
 
double dirichlet_gradient (Point point, scalar s, scalar cs,
         coord n, coord p, double bc, double * coef)
{
#if dimension == 2
  for (int _d = 0; _d < dimension; _d++)
    if (fabs(n.x) >= fabs(n.y))
      return dirichlet_gradient_x (point, s, cs, n, p, bc, coef);
#else // dimension == 3
  if (fabs(n.x) >= fabs(n.y)) {
    if (fabs(n.x) >= fabs(n.z))
      return dirichlet_gradient_x (point, s, cs, n, p, bc, coef);
  }
  else if (fabs(n.y) >= fabs(n.z))
    return dirichlet_gradient_y (point, s, cs, n, p, bc, coef);
  return dirichlet_gradient_z (point, s, cs, n, p, bc, coef);
#endif // dimension == 3
  return nodata;
}
 
bid embed;
 
/**
## Surface force and vorticity
 
We first define a function which computes
\f$\mathbf{\nabla}\mathbf{u}\cdot\mathbf{n}\f$ while taking the boundary
conditions on the embedded surface into account. */
 
static inline
coord embed_gradient (Point point, vector u, coord p, coord n)
{
  coord dudn;
  for (int _d = 0; _d < dimension; _d++) {
    bool dirichlet = false;
    double vb = u.x.boundary[embed] (point, point, u.x, &dirichlet);
    if (dirichlet) {
      double val;
      dudn.x = dirichlet_gradient (point, u.x, cs, n, p, vb, &val);
    }
    else // Neumann
      dudn.x = vb;
    if (dudn.x == nodata)
      dudn.x = 0.;
  }
  return dudn;
}
 
/**
The force exerted by the fluid on the solid can be written
\f[
\mathbf{F}_{\Gamma} = - \int_{\partial \Gamma} ( - p\mathbf{I} +
2 \mu \mathbf{D}) \cdot \mathbf{n}d \partial \Gamma
\f]
with \f$\Gamma\f$ the solid boundary. It can be further decomposed into a
pressure (i.e. "form") drag
\f[
\mathbf{F}_p = \int_{\partial \Gamma} p \mathbf{n}d \partial \Gamma
\f]
and a viscous drag
\f[
\mathbf{F}_{\mu} = - \int_{\partial \Gamma} 
2 \mu \mathbf{D} \cdot \mathbf{n}d \partial \Gamma
\f]
These two vectors are computed by the *embed_force()* function.
*/
 
trace
void embed_force (scalar p, vector u, vector mu, coord * Fp, coord * Fmu)
{
  coord Fps = {0}, Fmus = {0};
  foreach (reduction(+:Fps) reduction(+:Fmus), nowarning)
    if (cs[] > 0. && cs[] < 1.) {
 
      /**
      To compute the pressure force, we first get the coordinates of
      the barycentre of the embedded fragment, its area and normal,
      and then interpolate the pressure field on the surface. */
      
      coord n, b;
      double area = embed_geometry (point, &b, &n);
      area *= pow (Delta, dimension - 1);
      double Fn = area*embed_interpolate (point, p, b);
      for (int _d = 0; _d < dimension; _d++)
  Fps.x += Fn*n.x;
 
      /**
      To compute the viscous force, we first need to retrieve the
      local value of the viscosity (ideally at the barycentre of the
      embedded fragment). This is not completely trivial since it is
      defined on the faces of the cell. We use a
      surface-fraction-weighted average value. */
      
      if (constant(mu.x) != 0.) {
  double mua = 0., fa = 0.;
  for (int _d = 0; _d < dimension; _d++) {
    mua += mu.x[] + mu.x[1];
    fa  += fm.x[] + fm.x[1];
  }
  mua /= fa;
 
  /**
  To compute the viscous force, we need to take into account the
  (Dirichlet) boundary conditions for the velocity on the
  surface. We only know how to do this when computing the normal
  gradient \f$\mathbf{\nabla}\mathbf{u}\cdot\mathbf{n}\f$ using the
  [embed_gradient()](#embed_gradient) function. We thus
  need to re-express the viscous force using only normal
  derivatives of the velocity field.
  
  If we assume that \f$\mathbf{u}\f$ is constant on the boundary, then
  \f[
  \mathbf{{\nabla}} \mathbf{u} \cdot \mathbf{t}= \mathbf{0}
  \f]
  with \f$\mathbf{t}\f$ the unit tangent vector to the boundary. We
  thus have the relations
  \f[
  \mathbf{{\nabla}} \mathbf{u} = \left( \mathbf{{\nabla}} \mathbf{u}
  \cdot \mathbf{n} \right) \mathbf{n} + \left( \mathbf{{\nabla}}
  \mathbf{u} \cdot \mathbf{t} \right) \mathbf{t} = \left(
  \mathbf{{\nabla}} \mathbf{u} \cdot \mathbf{n} \right) \mathbf{n}
  \f]
  \f[
  \mathbf{D}= \frac{1}{2}  \left( \mathbf{{\nabla}} \mathbf{u} +
  \mathbf{{\nabla}}^T \mathbf{u} \right) = \frac{1}{2} 
  \left(\begin{array}{cc}
  2 \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_x & \left(
  \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_y + \left(
  \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x\\
  \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_y + \left(
  \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x & 2 \left(
  \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_y
  \end{array}\right)
  \f]
  \f[
  \mathbf{F}_{\mu} = - \int_{\Gamma} \left(\begin{array}{c}
  \left[2 \mu \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right)
  n_x \right] n_x + \mu \left[ \left( \mathbf{{\nabla}} u \cdot \mathbf{n}
  \right) n_y + \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x
  \right] n_y\\
  \left[2 \mu \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right)
  n_y \right] n_y + \mu \left[ \left( \mathbf{{\nabla}} u \cdot \mathbf{n}
  \right) n_y + \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x
  \right] n_x
  \end{array}\right)
  \f]
  \f[
  \mathbf{F}_{\mu} = - \int_{\Gamma} \left(\begin{array}{c}
  \mu \left[ \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) 
  (n^2_x + 1) + \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x
  n_y \right]\\
  \mu \left[ \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) 
  (n^2_y + 1) + \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_x
  n_y \right]
  \end{array}\right)
  \f]
  */
 
  assert (dimension == 2);
  coord dudn = embed_gradient (point, u, b, n);
  for (int _d = 0; _d < dimension; _d++)
    Fmus.x -= area*mua*(dudn.x*(sq(n.x) + 1.) + dudn.y*n.x*n.y);
      }
    }
 
  *Fp = Fps; *Fmu = Fmus;
}
 
/**
In two dimensions, *embed_vorticity()* returns the vorticity of
velocity field *u*, on the surface of the embedded boundary contained
in the cell. *p* is the relative position of the barycentre of the
embedded fragment and *n* its normal. */
 
#if dimension == 2
double embed_vorticity (Point point, vector u, coord p, coord n)
{
  /**
  We compute \f$\mathbf{{\nabla}}\mathbf{u}\cdot\mathbf{n}\f$, taking
  the boundary conditions into account. */
    
  coord dudn = embed_gradient (point, u, p, n);
 
  /**
  The vorticity is then obtained using the relations
  \f[
  \omega = \partial_x v - \partial_y u = 
  \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x - 
  \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_y
  \f]
  */
    
  return dudn.y*n.x - dudn.x*n.y;
}
#endif // dimension == 2
 
/**
## Flux through the embedded boundary
 
This function computes the flux through the embedded boundary contained 
within a cell
\f[
\int_b \mu \nabla s\cdot\mathbf{n} db
\f]
with \f$db\f$ the elementary boundary surface and \f$\mathbf{n}\f$ the embedded
boundary (outward-pointing) normal.
 
Boundary conditions for *s* are taken into account.
 
The result is returned in *val*. 
 
For degenerate cases, the value returned by the function must be
multiplied by `s[]` and added to *val*. */
 
double embed_flux (Point point, scalar s, vector mu, double * val)
{
 
  /**
  If the cell does not contain a fragment of embedded boundary, the
  flux is zero. */
  
  *val = 0.;
  if (cs[] >= 1. || cs[] <= 0.)
    return 0.;
 
  /**
  If the boundary condition is homogeneous Neumann, the flux is
  zero. */
 
  bool dirichlet = false;
  double grad = s.boundary[embed] (point, point, s, &dirichlet);
  if (!grad && !dirichlet)
    return 0.;
 
  /**
  We compute the normal, area and barycenter of the fragment of embedded
  boundary contained within the cell. */
 
  coord n = facet_normal (point, cs, fs), p;
  double alpha = plane_alpha (cs[], n);
  double area = plane_area_center (n, alpha, &p);
  if (metric_embed_factor)
    area *= metric_embed_factor (point, p);
 
  /**
  If the boundary condition is Dirichlet, we need to compute the
  normal gradient. */
 
  double coef = 0.;
  if (dirichlet) {
    normalize (&n);
    grad = dirichlet_gradient (point, s, cs, n, p, grad, &coef);
  }
 
  /**
  We retrieve the (average) value of \f$\mu\f$ without the metric. */
  
  double mua = 0., fa = 0.;
  for (int _d = 0; _d < dimension; _d++) {
    mua += mu.x[] + mu.x[1];
    fa  += fm.x[] + fm.x[1];
  }
  *val = - mua/(fa + SEPS)*grad*area/Delta;
  return - mua/(fa + SEPS)*coef*area/Delta;
}
 
/**
For ease of use, we replace the Neumann and Dirichlet functions with
macros so that they can be used either for standard domain boundaries
or for embedded boundaries. The distinction between the two cases is
based on whether the `dirichlet` parameter is passed to the boundary
function (using the `data` parameter). */
 
macro2
double dirichlet (double expr, Point point = point,
      scalar s = _s, bool * data = data)
{
  return data ? embed_area_center (point, &x, &y, &z),
    *((bool *)data) = true, expr : 2.*expr - s[];
}
 
macro2
double dirichlet_homogeneous (double expr, Point point = point,
            scalar s = _s, bool * data = data)
{
  return data ? *((bool *)data) = true, 0 : - s[];
}
 
macro2
double neumann (double expr, Point point = point,
    scalar s = _s, bool * data = data)
{
  return data ? embed_area_center (point, &x, &y, &z),
    *((bool *)data) = false, expr : Delta*expr + s[];
}
 
macro2
double neumann_homogeneous (double expr, Point point = point,
          scalar s = _s, bool * data = data)
{
  return data ? *((bool *)data) = false, 0 : s[];
}
 
/**
## Prolongation for the multigrid solver
 
We use a simplified prolongation operator for the [multigrid
solver](poisson.h#mg_solve) i.e. simple injection if bilinear
interpolation would use values which are fully contained within the
embedded boundary. */
 
#if MULTIGRID
static inline double bilinear_embed (Point point, scalar s)
{
  if (!coarse(cs) || !coarse(cs,child.x))
    return coarse(s);
  #if dimension >= 2
  if (!coarse(cs,0,child.y) || !coarse(cs,child.x,child.y))
    return coarse(s);
  #endif
  #if dimension >= 3
  if (!coarse(cs,0,0,child.z) || !coarse(cs,child.x,0,child.z) ||
      !coarse(cs,0,child.y,child.z) ||
      !coarse(cs,child.x,child.y,child.z))
    return coarse(s);  
  #endif
  return bilinear (point, s);
}
 
#define bilinear(point, s) bilinear_embed(point, s)
#endif // MULTIGRID
 
/**
## Lifting the "small cell" CFL restriction
 
For explicit advection schemes, the timestep is limited by the CFL
conditions 
\f[ 
\Delta t < \frac{c_s\Delta}{f_i|u_i|}
\f] 
where \f$i\f$ is the index of each face, and \f$c_s\f$ and \f$f_i\f$ are the
embedded volume and face fractions respectively. It is clear that the
timestep may need to be arbitrarily small if \f$c_s/f_s\f$ tends toward
zero. This is the "small cell" restriction of cut-cell finite-volume
techniques.
 
A classical technique to avoid this limitation is to use a "cell
merging" procedure, where the fluxes from cells which would "overflow"
are redistributed to neighboring cells.
 
The function below uses this approach to update a field *f*, advected
by the face velocity field *uf*, with corresponding advection fluxes
*flux*, during timestep *dt* which only verifies the standard CFL
condition
\f[
\Delta t < \frac{\Delta}{|u_i|}
\f]
*/
 
trace
void update_tracer (scalar f, vector uf, vector flux, double dt)
{
 
  /**
  Note that the distinction should be made between \f$c_m\f$, the cell
  fraction metric, and \f$c_s\f$, the embedded fraction. This is not done
  now so that embedded boundaries cannot be combined with a metric
  yet. 
 
  The field *e* will store the "overflowing" sum of fluxes for each cell. */
  
  scalar e[];
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
 
    /**
    If the cell is empty, it cannot overflow. */
    
    if (cs[] <= 0.)
      e[] = 0.;
 
    /**
    If the cell does not contain an embedded boundary, it cannot
    overflow either and the sum of the fluxes can be added to advance
    *f* in time. */
    
    else if (cs[] >= 1.) {
      for (int _d = 0; _d < dimension; _d++)
  f[] += dt*(flux.x[] - flux.x[1])/Delta;
      e[] = 0.;
    }
 
    /**
    If the cell contains the embedded boundary, we compute the maximum
    timestep verifying the restrictive CFL condition
    \f[ 
    \Delta t_{max} = \frac{c_s\Delta}{max(f_i|u_i|)}
    \f]
    Note that *fs* does not appear in the code below because *uf*
    already stores the product \f$u_f \, f_m\f$. */
    
    else {
      double umax = 0.;
      for (int i = 0; i <= 1; i++)
  for (int _d = 0; _d < dimension; _d++)
    if (fabs(uf.x[i]) > umax)
      umax = fabs(uf.x[i]);
      double dtmax = Delta*cm[]/(umax + SEPS);
 
      /**
      We compute the sum of the fluxes. */
      
      double F = 0.;
      for (int _d = 0; _d < dimension; _d++)
  F += flux.x[] - flux.x[1];
      F /= Delta*cm[];
 
      /**
      If the timestep is smaller than \f$\Delta t_{max}\f$, the cell
      cannot overflow and *f* is advanced in time using the entire
      flux. */
      
      if (dt <= dtmax) {
  f[] += dt*F;
  e[] = 0.;
      }
      
      /**
      Otherwise, the cell is filled "to the brim" by advancing *f*
      using the maximum allowable timestep. The field *e* is used to
      store the excess flux, weighted by the sum of the neighboring
      embedded fractions. */
      
      else {
  f[] += dtmax*F;
  double scs = 0.;
  for (int _n = 0; _n < 1; _n++)
    scs += sq(cm[]);  //Correct?
  e[] = (dt - dtmax)*F*cm[]/scs;
      }
    }
  }
 
  /**
  In a second phase, the excesses in each cell are added to the
  neighboring cells in proportion of their embedded fractions. */ 
 
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    double se = 0.;
    for (int _n = 0; _n < 1; _n++)
      se += e[];
    f[] += cs[]*se;
  }
}
 
/**
## Default settings
 
To apply the volume/area fraction-weighting to the solvers, we define
the embedded solid using the embedded fractions. The fields are no
longer constant and must be allocated.
 
The *cm* and *fm* fields contain the product of the metric and solid
factors. For a Cartesian coordinate system *cm* and *fm* are thus
identical to the solid factor fields *cs* and *fs*. */
 
/** @brief Event: metric (i = 0) */
void event_metric (void) 
{
  if (is_constant (fm.x)) {
    for (int _d = 0; _d < dimension; _d++)
      assert (constant (fm.x) == 1.);
    fm = fs;
  }
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    fs.x[] = 1.;
  if (is_constant (cm)) {
    assert (constant (cm) == 1.);
    cm = cs;
  }
  for (int _i = 0; _i < _N; _i++) /* foreach */
    cs[] = 1.;
 
#if TREE
  cs.refine = embed_fraction_refine;
 
  /**
  For prolongation we cannot use the same function since the surface
  fraction field *fs* is not necessarily defined for prolongation
  cells. So we switch back to the default fraction refinement (which
  is less accurate but only relies on *cs*). */
 
  // fixme: could work better with embed_fraction_refine
  // see porous1.tst
  cs.prolongation = fraction_refine;
  for (int _d = 0; _d < dimension; _d++)
    fs.x.prolongation = embed_face_fraction_refine_x;
  
  /**
  Note that we do not need to change the `refine` method since the
  default `refine` method calls the prolongation method for each
  component. */
  
#endif
  restriction ({cs, fs});
}
 
/**
We add the embedded boundary to the default display. */
 
/** @brief Event: defaults (i = 0 ) */
void event_defaults (void) {
  display ("draw_vof (c = 'cs', s = 'fs', filled = -1, "
     "fc = {0.5,0.5,0.5}, order = 2);");
}
 
/**
## References
 
~~~bib
@article{johansen1998,
  title={A Cartesian grid embedded boundary method for Poisson's
  equation on irregular domains},
  author={Johansen, Hans and Colella, Phillip},
  journal={Journal of Computational Physics},
  volume={147},
  number={1},
  pages={60--85},
  year={1998},
  publisher={Elsevier},
  url={https://pdfs.semanticscholar.org/17cd/babecd054d58da05c2ba009cccb3c687f58f.pdf}
}
 
@article{schwartz2006,
  title={A Cartesian grid embedded boundary method for the heat equation 
  and Poisson’s equation in three dimensions},
  author={Schwartz, Peter and Barad, Michael and Colella, Phillip and Ligocki, 
  Terry},
  journal={Journal of Computational Physics},
  volume={211},
  number={2},
  pages={531--550},
  year={2006},
  publisher={Elsevier},
  url={https://cloudfront.escholarship.org/dist/prd/content/qt0fp606kk/qt0fp606kk.pdf}
}
~~~
 
## See also
 
* [Notes on drag force computation](/src/notes/drag.tm)
*/