Interpolation

By default, probes and profile values are "P0" interpolated, that is their value is that of the containing cell or face, or closest vertex.

For cell-based values, a "P1" piecewise linear interpolation may be defined, using gradient reconstruction..

The following example shows how intepolation may be used in the cs_user_postprocess_values function.

First, we define an interpolation function:

 
/*----------------------------------------------------------------------------
 * Interpolate values defined on a mesh location at a given set of
 * points using a P0 interpolation.
 *
 * This function assumes the input is a field name. If no field matches
 * the name, a "generic" interpolation (assuming homogeneous Neumann boundary
 * conditions) is used).
 *
 * \param[in, out]  input           pointer to optional (untyped) value
 *                                  or structure.
 * \param[in]       datatype        associated datatype
 * \param[in]       val_dim         dimension of data values
 * \param[in]       n_points        number of interpolation points
 * \param[in]       point_location  location of points in mesh elements
 * \param[in]       point_coords    point coordinates
 * \param[in]       location_vals   values at mesh location
 * \param[out]      point_vals      interpolated values at points
 *----------------------------------------------------------------------------*/
 
static void
_cs_interpolate_from_location_p1(void                *input,
                                 cs_datatype_t        datatype,
                                 int                  val_dim,
                                 cs_lnum_t            n_points,
                                 const cs_lnum_t      point_location[],
                                 const cs_real_3_t    point_coords[],
                                 const void          *location_vals,
                                 void                *point_vals)
{
  /* If used with a non-real argument type, use P0 interpolation */
 
  if (   datatype != CS_REAL_TYPE
      || (val_dim != 1 && val_dim != 3 && val_dim != 6)) {
    cs_interpolate_from_location_p0(input,
                                    datatype,
                                    val_dim,
                                    n_points,
                                    point_location,
                                    point_coords,
                                    location_vals,
                                    point_vals);
    return;
  }
 
  /* Main usage */
 
  cs_mesh_t *m = cs_glob_mesh;
  const cs_real_3_t *cell_cen
    = (const cs_real_3_t *)cs_glob_mesh_quantities->cell_cen;
 
  cs_gradient_type_t gradient_type = CS_GRADIENT_ITER;
  cs_halo_type_t halo_type = CS_HALO_STANDARD;
 
  cs_gradient_type_by_imrgra(cs_glob_space_disc->imrgra,
                             &gradient_type,
                             &halo_type);
 
  cs_field_t *f = NULL;
  if (input != NULL) {
    const char *name = input;
    f = cs_field_by_name_try(name);
  }
 
  switch(val_dim) {
  case 1:
    {
      const cs_real_t *c_vals = (const cs_real_t *)location_vals;
      cs_real_t *p_vals = (cs_real_t *)point_vals;
 
      cs_real_3_t *grad;
      BFT_MALLOC(grad, m->n_cells_with_ghosts, cs_real_3_t);
      if (f != NULL)
        cs_field_gradient_scalar(f,
                                 false, /* use_previous_t */
                                 1, /* inc */
                                 true, /* _recompute_cocg */
                                 grad);
      else
        cs_gradient_scalar_synced_input("user",
                                        gradient_type,
                                        halo_type,
                                        1,      /* inc */
                                        true,   /* recompute_cocg */
                                        100,    /* n_r_sweeps */
                                        0,      /* tr_dim */
                                        0,      /* hyd_p_flag */
                                        0,      /* w_stride */
                                        0,      /* verbosity */
                                        -1,     /* clip_mode */
                                        1e-5,   /* epsilon */
                                        0,      /* extrap */
                                        1.5,    /* clip_coeff */
                                        NULL,   /* f_ext[] */
                                        NULL,   /* bc_coeff_a */
                                        NULL,   /* bc_coeff_b */
                                        c_vals,
                                        NULL,   /* c_weight */
                                        NULL,   /* cpl */
                                        grad);
 
      for (cs_lnum_t i = 0; i < n_points; i++) {
        cs_lnum_t c_id = point_location[i];
        if (c_id > -1) {
          cs_real_t d[3] = {point_coords[i][0] - cell_cen[c_id][0],
                            point_coords[i][1] - cell_cen[c_id][1],
                            point_coords[i][2] - cell_cen[c_id][2]};
 
          p_vals[i] = c_vals[c_id] + grad[c_id][0]*d[0]
                                   + grad[c_id][1]*d[1]
                                   + grad[c_id][2]*d[2];
        }
        else
          p_vals[i] = 0;
      }
 
      BFT_FREE(grad);
    }
    break;
 
  case 3:
    {
      const cs_real_3_t *c_vals = (const cs_real_3_t *)location_vals;
      cs_real_3_t *p_vals = (cs_real_3_t *)point_vals;
 
      cs_real_33_t *grad;
      BFT_MALLOC(grad, m->n_cells_with_ghosts, cs_real_33_t);
      if (f != NULL)
        cs_field_gradient_vector(f,
                                 false, /* use_previous_t */
                                 1, /* inc */
                                 grad);
      else
        cs_gradient_vector_synced_input("user",
                                        gradient_type,
                                        halo_type,
                                        1,      /* inc */
                                        100,    /* n_r_sweeps */
                                        0,      /* verbosity */
                                        -1,     /* clip_mode */
                                        1e-5,   /* epsilon */
                                        1.5,    /* clip_coeff */
                                        NULL,   /* bc_coeff_a */
                                        NULL,   /* bc_coeff_b */
                                        c_vals,
                                        NULL,   /* c_weight */
                                        NULL,   /* cpl */
                                        grad);
 
      for (cs_lnum_t i = 0; i < n_points; i++) {
        cs_lnum_t c_id = point_location[i];
        if (c_id > -1) {
          cs_real_t d[3] = {point_coords[i][0] - cell_cen[c_id][0],
                            point_coords[i][1] - cell_cen[c_id][1],
                            point_coords[i][2] - cell_cen[c_id][2]};
 
          for (cs_lnum_t j = 0; j < 3; j++) {
            p_vals[i][j] = c_vals[c_id][j] + grad[c_id][j][0]*d[0]
                                           + grad[c_id][j][1]*d[1]
                                           + grad[c_id][j][2]*d[2];
          }
        }
        else {
          for (cs_lnum_t j = 0; j < 6; j++)
            p_vals[i][j] = 0;
        }
      }
 
      BFT_FREE(grad);
    }
    break;
 
  case 6:
    {
      const cs_real_6_t *c_vals = (const cs_real_6_t *)location_vals;
      cs_real_6_t *p_vals = (cs_real_6_t *)point_vals;
 
      cs_real_63_t *grad;
      BFT_MALLOC(grad, m->n_cells_with_ghosts, cs_real_63_t);
      if (f != NULL)
        cs_field_gradient_tensor(f,
                                 false, /* use_previous_t */
                                 1, /* inc */
                                 grad);
      else
        cs_gradient_tensor_synced_input("user",
                                        gradient_type,
                                        halo_type,
                                        1,      /* inc */
                                        100,    /* n_r_sweeps */
                                        0,      /* verbosity */
                                        -1,     /* clip_mode */
                                        1e-5,   /* epsilon */
                                        1.5,    /* clip_coeff */
                                        NULL,   /* bc_coeff_a */
                                        NULL,   /* bc_coeff_b */
                                        c_vals,
                                        grad);
 
      for (cs_lnum_t i = 0; i < n_points; i++) {
        cs_lnum_t c_id = point_location[i];
        if (c_id > -1) {
          cs_real_t d[3] = {point_coords[i][0] - cell_cen[c_id][0],
                            point_coords[i][1] - cell_cen[c_id][1],
                            point_coords[i][2] - cell_cen[c_id][2]};
 
          for (cs_lnum_t j = 0; j < 6; j++) {
            p_vals[i][j] = c_vals[c_id][j] + grad[c_id][j][0]*d[0]
                                           + grad[c_id][j][1]*d[1]
                                           + grad[c_id][j][2]*d[2];
          }
        }
        else {
          for (cs_lnum_t j = 0; j < 6; j++)
            p_vals[i][j] = 0;
        }
      }
 
      BFT_FREE(grad);
    }
    break;
 
  default:
    assert(0);
  }
}
 

Note the the gradient reconstruction used here assumes ghost cell values are synchronized, which should be the case for field values at this calling stage.

This interpolation function may the be passed to the probe values output function (in cs_user_postprocess_values):

    const char *name = cs_probe_set_get_name(probes);
 
    int n_p_fields = 2;
    const char *p_field_names[] = {"velocity", "temperature"};
 
    for (int i = 0; i < n_p_fields; i++) {
 
      cs_field_t *f = cs_field_by_name_try(p_field_names[i]);
 
      if (f != NULL) {
 
        /* use different name to avoid conflict with field name in case already
           present in probe set through default output */
 
        char p_name[64];
        snprintf(p_name, 63, "%s_p", f->name); p_name[63] = '\0';
 
        cs_post_write_probe_values
          (mesh_id,
           CS_POST_WRITER_ALL_ASSOCIATED,    /* writer id filter */
           p_name,                           /* var_name */
           f->dim,                           /* var_dim */
           CS_POST_TYPE_cs_real_t,
           1,                                /* parent location id */
           _cs_interpolate_from_location_p1, /* P1 interpolation */
           f->name,                          /* interpolation input */
           f->val,
           ts);
 
      }
 
    }
 
  }

In this case, selected outputs are named by appending "_p" to the field name to allow combining default probe outputs with interpolated outputs of specific fields.

For simplicity here, values are output to the main probe set and writer, which is assumed to be defined using the GUI in this example.

Note also that interpolation could be also used in some cs_user_extra_operations cases.