Advanced definitions of post-processing and mesh zones

Output meshes defined by user-defined functions

More advanced mesh element selection is possible using cs_post_define_volume_mesh_by_func or cs_post_define_surface_mesh_by_func, which allow defining volume or surface meshes using user-defined element lists.

The possibility to modify a mesh over time is limited by the most restrictive writer which is associated with. For instance, if writer 1 allows the modification of the mesh topology (argument time_dep = FVM_WRITER_TRANSIENT_CONNECT in the call to cs_post_define_writer) and writer 2 allows no modification (time_dep = FVM_WRITER_FIXED_MESH), a user post-processing mesh associated with writers 1 and 2 will not be modifiable, but a mesh associated with writer 1 only will be modifiable. The modification can be done by using the advanced cs_post_define_volume_mesh_by_func or cs_post_define_surface_mesh_by_func, associated with a user-defined selection function based on time-varying criteria (such as field values being above a given threshold). If the time_dep argument is set to true, the mesh will be redefined using the selection function at each output time step for every modifiable mesh.

Example: surface mesh with complex selection criteria

In the following example, we build a surface mesh containing interior faces separating cells of group "2" from those of group "3", (assuming no cell has both colors), as well as boundary faces of group "4".

This is done by first defining 2 selection functions, whose arguments and behavior match the cs_post_elt_select_t type.

The function for selection of interior faces separating cells of two groups also illustrates the usage of the cs_selector_get_family_list function to build a mask allowing direct checking of this criterion when comparing cells adjacent to a given face:

static void
_i_faces_select_example(void         *input,
                        cs_lnum_t    *n_faces,
                        cs_lnum_t   **face_ids)
{
  CS_UNUSED(input);
 
  cs_lnum_t i, face_id;
  cs_lnum_t n_families = 0;
  cs_int_t *family_list = NULL;
  int *family_mask = NULL;
 
  cs_lnum_t n_i_faces = 0;
  cs_lnum_t *i_face_ids = NULL;
 
  const cs_mesh_t *m = cs_glob_mesh;
 
  /* Allocate selection list */
 
  BFT_MALLOC(i_face_ids, m->n_i_faces, cs_lnum_t);
 
  /* Build mask on families matching groups "2" (1), "3" (2) */
 
  BFT_MALLOC(family_list, m->n_families, cs_int_t);
  BFT_MALLOC(family_mask, m->n_families, int);
 
  for (i = 0; i < m->n_families; i++)
    family_mask[i] = 0;
 
  cs_selector_get_family_list("2",  &n_families, family_list);
 
  for (i = 0; i < n_families; i++)
    family_mask[family_list[i] - 1] += 1;
 
  cs_selector_get_family_list("3",  &n_families, family_list);
 
  for (i = 0; i < n_families; i++)
    family_mask[family_list[i] - 1] += 2;
 
  BFT_FREE(family_list);
 
  /* Now that mask is built, test for adjacency */
 
  for (face_id = 0; face_id < m->n_i_faces; face_id++) {
 
    /* Adjacent cells  and flags */
 
    cs_lnum_t c1 = m->i_face_cells[face_id][0];
    cs_lnum_t c2 = m->i_face_cells[face_id][1];
 
    int iflag1 = family_mask[m->cell_family[c1]];
    int iflag2 = family_mask[m->cell_family[c2]];
 
    /* Should the face belong to the extracted mesh ? */
 
    if ((iflag1 == 1 && iflag2 == 2) || (iflag1 == 2 && iflag2 == 1)) {
      i_face_ids[n_i_faces] = face_id;
      n_i_faces += 1;
    }
 
  }
 
  /* Free memory */
 
  BFT_FREE(family_mask);
  BFT_REALLOC(i_face_ids, n_i_faces, cs_lnum_t);
 
  /* Set return values */
 
  *n_faces = n_i_faces;
  *face_ids = i_face_ids;
}

The function for selection of boundary faces is simpler, as it simply needs to apply the selection criterion for boundary faces:

static void
_i_faces_select_example(void         *input,
                        cs_lnum_t    *n_faces,
                        cs_lnum_t   **face_ids)
{
  CS_UNUSED(input);
 
  cs_lnum_t i, face_id;
  cs_lnum_t n_families = 0;
  cs_int_t *family_list = NULL;
  int *family_mask = NULL;
 
  cs_lnum_t n_i_faces = 0;
  cs_lnum_t *i_face_ids = NULL;
 
  const cs_mesh_t *m = cs_glob_mesh;
 
  /* Allocate selection list */
 
  BFT_MALLOC(i_face_ids, m->n_i_faces, cs_lnum_t);
 
  /* Build mask on families matching groups "2" (1), "3" (2) */
 
  BFT_MALLOC(family_list, m->n_families, cs_int_t);
  BFT_MALLOC(family_mask, m->n_families, int);
 
  for (i = 0; i < m->n_families; i++)
    family_mask[i] = 0;
 
  cs_selector_get_family_list("2",  &n_families, family_list);
 
  for (i = 0; i < n_families; i++)
    family_mask[family_list[i] - 1] += 1;
 
  cs_selector_get_family_list("3",  &n_families, family_list);
 
  for (i = 0; i < n_families; i++)
    family_mask[family_list[i] - 1] += 2;
 
  BFT_FREE(family_list);
 
  /* Now that mask is built, test for adjacency */
 
  for (face_id = 0; face_id < m->n_i_faces; face_id++) {
 
    /* Adjacent cells  and flags */
 
    cs_lnum_t c1 = m->i_face_cells[face_id][0];
    cs_lnum_t c2 = m->i_face_cells[face_id][1];
 
    int iflag1 = family_mask[m->cell_family[c1]];
    int iflag2 = family_mask[m->cell_family[c2]];
 
    /* Should the face belong to the extracted mesh ? */
 
    if ((iflag1 == 1 && iflag2 == 2) || (iflag1 == 2 && iflag2 == 1)) {
      i_face_ids[n_i_faces] = face_id;
      n_i_faces += 1;
    }
 
  }
 
  /* Free memory */
 
  BFT_FREE(family_mask);
  BFT_REALLOC(i_face_ids, n_i_faces, cs_lnum_t);
 
  /* Set return values */
 
  *n_faces = n_i_faces;
  *face_ids = i_face_ids;
}

Given these tow functions, the mesh can be defined using the cs_post_define_surface_mesh_by_func function, passing it the user-defined selection functions (actually, function pointers):

  {
    const int n_writers = 1;
    const int writer_ids[] = {1};  /* Associate to writer 1 */
 
    /* Define postprocessing mesh */
 
    cs_post_define_surface_mesh_by_func(3,               /* mesh id */
                                        "Mixed surface",
                                        _i_faces_select_example,
                                        _b_faces_select_example,
                                        NULL,            /* i_faces_sel_input */
                                        NULL,            /* b_faces_sel_input */
                                        false,           /* time varying */
                                        false,           /* add_groups */
                                        false,           /* auto_variables */
                                        n_writers,
                                        writer_ids);
  }

Example: time-varying mesh

A mesh defined through the advanced cs_post_define_surface_mesh_by_func, cs_post_define_volume_mesh_by_func, or cs_post_define_particles_mesh_by_func may vary in time, as long as the matching time_varying argument is set to true, and the mesh (or aliases thereof) id only associated to writers defined with the FVM_WRITER_TRANSIENT_CONNECT option. In the case of particles, which always vary in time, this allows also varying the selection (filter) function with time.

In the following example, we build a volume mesh containing cells with values of field named "He_fraction" greater than 0.05.

First, we define the selection function:

static void
_he_fraction_05_select(void        *input,
                       cs_lnum_t   *n_cells,
                       cs_lnum_t  **cell_ids)
{
  CS_UNUSED(input);
 
  cs_lnum_t _n_cells = 0;
  cs_lnum_t *_cell_ids = NULL;
 
  const cs_mesh_t *m = cs_glob_mesh;
 
  cs_field_t *f = cs_field_by_name("He_fraction"); /* Get access to field */
 
  if (f == NULL)
    bft_error(__FILE__, __LINE__, 0,
              "No field with name \"He_fraction\" defined");
 
  /* Before time loop, field is defined, but has no values yet,
     so ignore that case (postprocessing mesh will be initially empty) */
 
  if (f->val != NULL) {
 
    BFT_MALLOC(_cell_ids, m->n_cells, cs_lnum_t); /* Allocate selection list */
 
    for (cs_lnum_t i = 0; i < m->n_cells; i++) {
      if (f->val[i] > 5.e-2) {
        _cell_ids[_n_cells] = i;
        _n_cells += 1;
      }
    }
 
    BFT_REALLOC(_cell_ids, _n_cells, cs_lnum_t); /* Adjust size (good practice,
                                                    but not required) */
 
  }
 
  /* Set return values */
 
  *n_cells = _n_cells;
  *cell_ids = _cell_ids;
}

Then, we simply define matching volume mesh passing the associated selection function pointer:

  {
    const int n_writers = 1;
    const int writer_ids[] = {2};  /* Associate to writer 2 */
 
    /* Define postprocessing mesh */
 
    cs_post_define_volume_mesh_by_func(4,               /* mesh id */
                                       "He_fraction_05",
                                       _he_fraction_05_select,
                                       NULL,            /* _c_05_select_input */
                                       true,            /* time varying */
                                       false,           /* add_groups */
                                       false,           /* auto_variables */
                                       n_writers,
                                       writer_ids);
  }

The matching function will be called at all time steps requiring output of this mesh.

Warning

some mesh formats do not allow changing meshes (or the implemented output functions do not allow them yet) and some may not allow empty meshes, even if this is only transient.

Other advanced mesh types

Example: edges mesh

In cases where a mesh containing polygonal elements is output through a writer configured to divide polygons into triangles (for example when visualization tools do not support polygons, or when highly non convex faces lead to visualization artifacts), it may be useful to extract a mesh containing the edges of the original mesh so as to view the polygon boundaries as an overlay.

In the following example, we build such a mesh (with id 5), based on the faces of a mesh with id 1:

  {
    const int n_writers = 1;
    const int writer_ids[] = {4};  /* Associate to writer 4 */
 
    cs_post_define_edges_mesh(5, /* mesh_id */
                              1, /* base_mesh_id */
                              n_writers,
                              writer_ids);
  }

Output of various space-filling curves

In the below example, edge meshes illustrating various space-filling curves possibilities are output.

First the below functions write the various space-filling curves either in serial or parallel.

 
/*----------------------------------------------------------------------------
 * Write space-filling curves for main mesh
 *
 * parameters:
 *   writer <-- FVM writer
 *---------------------------------------------------------------------------*/
 
static void
_cs_post_write_sfc_serial(fvm_writer_t   *writer)
{
  fvm_io_num_sfc_t  sfc_id;
  cs_lnum_t  i, j, k;
 
  cs_lnum_t *connect = NULL, *order = NULL;
  double *coords = NULL, *val = NULL;
  fvm_nodal_t *nm = NULL;
  fvm_io_num_t *io_num = NULL;
 
  const cs_mesh_t *m = cs_glob_mesh;
  const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
  const cs_lnum_t n_edges = m->n_cells - 1;
  const cs_gnum_t *cell_gnum = NULL;
  const double *var_ptr[1] = {NULL};
 
  BFT_MALLOC(order, m->n_cells, cs_lnum_t);
  BFT_MALLOC(val, m->n_cells, double);
  BFT_MALLOC(coords, m->n_cells*3, double);
 
  /* Loop on space-filling curve types */
 
  for (sfc_id = FVM_IO_NUM_SFC_MORTON_BOX;
       sfc_id <= FVM_IO_NUM_SFC_HILBERT_CUBE;
       sfc_id++) {
 
    BFT_MALLOC(connect, n_edges*2, cs_lnum_t);
 
    io_num = fvm_io_num_create_from_sfc(mq->cell_cen,
                                        3,
                                        m->n_cells,
                                        sfc_id);
 
    cell_gnum = fvm_io_num_get_global_num(io_num);
 
    cs_order_gnum_allocated(NULL, cell_gnum, order, m->n_cells);
 
    for (i = 0; i < m->n_cells; i++) {
      j = order[i];
      for (k = 0; k < 3; k++)
        coords[i*3 + k] = mq->cell_cen[j*3 + k];
      val[i] = i+1;
    }
 
    for (i = 0; i < n_edges; i++) {
      connect[i*2] = i+1;
      connect[i*2+1] = i+2;
    }
 
    cell_gnum = NULL;
    fvm_io_num_destroy(io_num);
 
    nm = fvm_nodal_create(fvm_io_num_sfc_type_name[sfc_id], 3);
 
    fvm_nodal_append_by_transfer(nm,
                                 m->n_cells - 1,
                                 FVM_EDGE,
                                 NULL,
                                 NULL,
                                 NULL,
                                 connect,
                                 NULL);
 
    fvm_nodal_set_shared_vertices(nm, coords);
 
    fvm_writer_export_nodal(writer, nm);
 
    var_ptr[0] = val;
 
    fvm_writer_export_field(writer,
                            nm,
                            _("order"),
                            FVM_WRITER_PER_NODE,
                            1,
                            CS_INTERLACE,
                            0,
                            0,
                            CS_DOUBLE,
                            -1,
                            0.0,
                            (const void * *)var_ptr);
 
    fvm_nodal_destroy(nm);
  }
 
  /* Free memory */
 
  BFT_FREE(order);
  BFT_FREE(val);
  BFT_FREE(coords);
}
 
#if defined(HAVE_MPI)
 
/*----------------------------------------------------------------------------
 * Write space-filling curves for main mesh
 *
 * parameters:
 *   writer <-- FVM writer
 *---------------------------------------------------------------------------*/
 
static void
_cs_post_write_sfc_parall(fvm_writer_t  *writer)
{
  fvm_io_num_sfc_t  sfc_id;
  cs_lnum_t  i;
 
  cs_lnum_t *connect = NULL, *order = NULL;
  cs_gnum_t *vtx_gnum = NULL, *edge_gnum = NULL;
  double *val = NULL;
  cs_coord_t *coords = NULL;
  fvm_nodal_t *nm = NULL;
  fvm_io_num_t *io_num = NULL;
 
  const cs_mesh_t *m = cs_glob_mesh;
  const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
  const cs_gnum_t *cell_gnum = NULL;
  const double  *var_ptr[1] = {NULL};
 
  /* Loop on space-filling curve types */
 
  for (sfc_id = FVM_IO_NUM_SFC_MORTON_BOX;
       sfc_id <= FVM_IO_NUM_SFC_HILBERT_CUBE;
       sfc_id++) {
 
    cs_lnum_t block_size = 0;
    cs_lnum_t n_edges = 0;
    cs_block_dist_info_t bi;
    cs_part_to_block_t *d = NULL;
 
    io_num = fvm_io_num_create_from_sfc(mq->cell_cen,
                                        3,
                                        m->n_cells,
                                        sfc_id);
 
    cell_gnum = fvm_io_num_get_global_num(io_num);
 
    /* Distribute to blocks so that edge connectivity is trivial */
 
    bi = cs_block_dist_compute_sizes(cs_glob_rank_id,
                                     cs_glob_n_ranks,
                                     0,
                                     0,
                                     m->n_g_cells);
 
    d = cs_part_to_block_create_by_gnum(cs_glob_mpi_comm,
                                        bi,
                                        m->n_cells,
                                        cell_gnum);
 
    block_size = (bi.gnum_range[1] - bi.gnum_range[0]);
 
    if (block_size > 0) {
      BFT_MALLOC(connect, block_size*2, cs_lnum_t);
      BFT_MALLOC(val, block_size+1, double);
      BFT_MALLOC(coords, (block_size+1)*3, double);
      BFT_MALLOC(vtx_gnum, block_size+1, cs_gnum_t);
      BFT_MALLOC(edge_gnum, block_size, cs_gnum_t);
    }
 
    /* Distribute blocks on ranks */
 
    cs_part_to_block_copy_array(d,
                                CS_DOUBLE,
                                3,
                                mq->cell_cen,
                                coords);
 
    cell_gnum = NULL;
    fvm_io_num_destroy(io_num);
 
    cs_part_to_block_destroy(&d);
 
    /* Add vertex for connectivity with next rank */
 
    if (block_size > 0) {
      MPI_Status status;
      int prev_rank = cs_glob_rank_id - bi.rank_step;
      int next_rank = cs_glob_rank_id + bi.rank_step;
      if (prev_rank < 0)
        prev_rank = MPI_PROC_NULL;
      if (bi.gnum_range[1] > m->n_g_cells)
        next_rank = MPI_PROC_NULL;
      MPI_Sendrecv(coords, 3, MPI_DOUBLE, prev_rank, 0,
                   coords + 3*block_size, 3, MPI_DOUBLE, next_rank, 0,
                   cs_glob_mpi_comm, &status);
    }
 
    for (i = 0; i < block_size; i++) {
      vtx_gnum[i] = bi.gnum_range[0] + i;
      if (vtx_gnum[i] < m->n_g_cells) {
        connect[n_edges*2] = i+1;
        connect[n_edges*2+1] = i+2;
        edge_gnum[n_edges] = vtx_gnum[i];
        n_edges++;
      }
      val[i] = vtx_gnum[i];
    }
    if (block_size > 0) {
      vtx_gnum[block_size] = bi.gnum_range[0] + block_size;
      val[block_size] = vtx_gnum[block_size];
    }
 
    BFT_FREE(order);
 
    nm = fvm_nodal_create(fvm_io_num_sfc_type_name[sfc_id], 3);
 
    fvm_nodal_append_by_transfer(nm,
                                 n_edges,
                                 FVM_EDGE,
                                 NULL,
                                 NULL,
                                 NULL,
                                 connect,
                                 NULL);
 
    connect = NULL;
 
    fvm_nodal_set_shared_vertices(nm, coords);
 
    fvm_nodal_init_io_num(nm, edge_gnum, 1);
    fvm_nodal_init_io_num(nm, vtx_gnum, 0);
 
    fvm_writer_export_nodal(writer, nm);
 
    var_ptr[0] = val;
 
    fvm_writer_export_field(writer,
                            nm,
                            _("order"),
                            FVM_WRITER_PER_NODE,
                            1,
                            CS_INTERLACE,
                            0,
                            0,
                            CS_DOUBLE,
                            -1,
                            0.0,
                            (const void * *)var_ptr);
 
    /* Free memory */
 
    if (block_size > 0) {
      BFT_FREE(val);
      BFT_FREE(coords);
      BFT_FREE(vtx_gnum);
      BFT_FREE(edge_gnum);
    }
 
    fvm_nodal_destroy(nm);
 
  }
}
 
#endif /* defined(HAVE_MPI) */
 

Then a fake cell selection function is used to call the writing of the space-filling curves at the correct step.

 
/*----------------------------------------------------------------------------
 * Cell selection function adapted to generate space-filling curves.
 *
 * This is a specific case, where we do not actually select any cells, but
 * generate a temporary, specific writer (based on writer -1 settings),
 * and edge meshes illustrating the various space-filling curve
 * possibilities are output using this writer.
 *
 * Doing this with this function allows executing these steps once the
 * mesh is preprocessed.
 *----------------------------------------------------------------------------*/
 
static void
_sfc_cell_select(void        *input,
                 cs_lnum_t   *n_cells,
                 cs_lnum_t  **cell_ids)
{
  *n_cells = 0;
  *cell_ids = NULL;
 
  fvm_writer_t *w = NULL;
 
  /* Create default writer */
 
  w = fvm_writer_init("SFC",
                      "postprocessing",
                      cs_post_get_default_format(),
                      cs_post_get_default_format_options(),
                      FVM_WRITER_FIXED_MESH);
 
#if defined(HAVE_MPI)
  if (cs_glob_n_ranks > 1)
    _cs_post_write_sfc_parall(w);
#endif
 
  if (cs_glob_n_ranks == 1)
    _cs_post_write_sfc_serial(w);
 
  fvm_writer_finalize(w);
}
 

Finally edge meshes are defined to illustrate the various possibilities.