radial.h

Go to the documentation of this file.

/** @file radial.h
 */
/**
# Radial/cylindrical coordinates
 
This implements the radial coordinate mapping illustrated below.
 
![Radial coordinate mapping](figures/radial.svg)
 
It works in 1D, 2D and 3D. The 3D version corresponds to cylindrical
coordinates since the \f$z\f$-coordinate is unchanged. 
 
The only parameter is \f$d\theta\f$, the total angle of the sector. */
 
double dtheta = pi/3.;
 
/**
For convenience we add definitions for the radial and angular
coordinates $(r, \theta)$. */
 
macro VARIABLES (Point point = point, int _ig = ig, int _jg = jg, int _kg = kg)
{
  VARIABLES (point, _ig, _jg, _kg);
  double r = x, theta = y*dtheta/L0;
  NOT_UNUSED(r); NOT_UNUSED(theta);
}
 
/** @brief Event: metric (i = 0) */
void event_metric (void) 
{
 
  /**
  We initialise the scale factors, taking care to first allocate the
  fields if they are still constant. */
 
  if (is_constant(cm)) {
    scalar * l = list_copy (all);
    cm = {0} /* new scalar */;
    free (all);
    all = list_concat ({cm}, l);
    free (l);
  }
  if (is_constant(fm.x)) {
    scalar * l = list_copy (all);
    fm = {0} /* new vector */;
    free (all);
    all = list_concat ((scalar *){fm}, l);
    free (l);
  }
 
  /**
  The area (in 2D) of a mapped element is the area of an annulus
  defined by the two radii \f$r-\Delta/2\f$ and \f$r+\Delta/2\f$, divided by
  the total number of sectors \f$N=2\pi L0/(d\theta\Delta)\f$, this gives
  \f[
  \frac{\pi [(r + \Delta / 2)^2 - (r - \Delta / 2)^2]}{2 \pi L 0 / (d \theta
  \Delta)} = \frac{2 \pi r \Delta}{2 \pi L 0 / (d \theta \Delta)} = \frac{rd
  \theta}{L 0} \Delta^2
  \f]
  By definition, the (area) metric factor `cm` is the mapped area divided
  by the unmapped area \f$\Delta^2\f$. */
 
  scalar cmv = cm;
  for (int _i = 0; _i < _N; _i++) /* foreach */
    cmv[] = r*dtheta/L0;
 
  /**
  It is important to set proper boundary conditions, in particular
  when refining the grid. */
  
  /* BC: cm[left] = dirichlet(r*dtheta/L0) */
  /* BC: cm[right] = dirichlet(r*dtheta/L0) */
  
  /**
  The (length) metric factor `fm` is the ratio of the mapped length of
  a face to its unmapped length \f$\Delta\f$. In the present case, it is
  unity for all dimensions except for the \f$x\f$ coordinates for which it
  is the ratio of the arclength by the unmapped length \f$\Delta\f$. We
  also set a small minimal value to avoid division by zero, in the
  case of a vanishing inner radius. */
  
  vector fmv = fm;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    fmv.x[] = 1.;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    fmv.x[] = max(r*dtheta/L0, 1e-20);
}