Setting porosity values from a CAD file

Using an appropriate CAD representation, it may be possible to compute the porosity fields using "common" or "cut" boolean operations between cells in the porous zones and the CAD volume. If the CAD shape represents the actual fluid volume, the common volume between the shape and the porous mesh section's cells will be used. If the CAD shape represents the solid volume, it should be cut from the pourous mesh section's cells.

The following example uses functions from the OpenCascade Technology (OCCT) libraries (https://dev.opencascade.org/doc/overview/html/index.html), and assumes they are installed on the user's machine. As these libraries form the backbone of the SALOME platform's CAD features, they are present with any generic or common SALOME install, and are available as a package in many Linux distributions.

Libraries to compile and link.

As the OCCT libraries provide a C++ API, interfacing them with code_saturne requires a file written in C++, so the following example uses a separate C++ file, cs_cad_intersect.cxx to to most of the work. This file is also a user example, and may be improved or modified.

To compile and link, additional compiler flags must be passed to code_saturne. They may be defined in <_ref cs_user_scripts.py, and the following values are recommended (adapting the paths to the local environment):

occ_include_path = "/opt/occ/7.4.0/include/opencascade"
occ_lib_paths = ("/opt/occ/7.4.0/lib", "/opt/gl2ps/1.4.0.1/lib")

occ_libs "-lTKMesh -lTKernel -lTKG2d -lTKG3d -lTKMath -lTKIGES  -lTKXSBase -lTKBin -lTKBool -lTKBO -lTKCDF -lTKBRep -lTKTopAlgo -lTKGeomAlgo -lTKGeomBase -lTKOffset -lTKPrim -lTKSTEP -lTKSTEPBase -lTKSTEPAttr -lTKHLR -lTKFeat"

domain.compile_cxxflags = "-std=c++11 -I" + occ_include_path;
domain.compile_libs = ""
for p in occ_lib_paths:
    domain.compile_libs += "-L" + p + " -Wl,-rpath -Wl," + p + " "
    domain.compile_libs += occ_libs

Function prototypes.

As the prototypes for the cs_cad_intersect.cxx are not part of the defaul code_saturne installation, do not forget to include the cs_cad_intersect.h header in the Local headers section of the cs_user_porosity.c file:

#include "cs_cad_intersect.h"

Zone selection

It is recommended to use the standard zone selection mechanism to select cells which may be intersected by the CAD shape, for example:

  const cs_zone_t  *z_poro = cs_volume_zone_by_name("Zone_1");

Local variables and initialization

A few local variables may allow a more concise syntax, and optional face porosity arrays may be declared as temporary work arrays if face factors are required::

  cs_mesh_t *m = cs_glob_mesh;
  cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;

  const cs_lnum_t n_i_faces = m->n_i_faces;
  const cs_lnum_t n_b_faces = m->n_b_faces;

  /* Temporary arrays for face porosity */

  cs_real_t  *cell_porosity = cs_field_by_name("porosity")->val;

  cs_real_t  *i_face_porosity, *b_face_porosity;

  BFT_MALLOC(i_face_porosity, m->n_i_faces, cs_real_t);
  BFT_MALLOC(b_face_porosity, m->n_b_faces, cs_real_t);

  for (cs_lnum_t i = 0; i < m->n_i_faces; i++)
    i_face_porosity[i] = 1;
  for (cs_lnum_t i = 0; i < m->n_b_faces; i++)
    b_face_porosity[i] = 1;

CAD intersection operation

The actual CAD intersection is done using the following call (adapting the file path to the actual file):

  cs_cad_intersect(m,
                   "Compound_1.step",
                   CS_CAD_INTERSECT_CUT,
                   z_poro->n_elts,
                   z_poro->elt_ids,
                   cell_porosity,
                   NULL,  /* modified cell centers */
                   i_face_porosity,
                   NULL,  /* modified interior face centers */
                   b_face_porosity,
                   NULL);  /* modified boundary face centers */

Computing porous quantities

Though the cell porosity is directly set by the call above, face quantities may be computed from the face porosities. The following code also handles the cleanup op local arrays.

  /* synchronize ghost cells */

  if (m->halo != NULL)
    cs_halo_sync_var(m->halo, CS_HALO_EXTENDED, cell_porosity);

  /* Set interior and boundary face values */

  if (mq->i_f_face_surf != mq->i_face_surf) {

    const cs_real_3_t *restrict i_face_normal
      = (const cs_real_3_t *restrict)mq->i_face_normal;
    cs_real_3_t *restrict i_f_face_normal
      = (cs_real_3_t *restrict)mq->i_f_face_normal;
    cs_real_t *i_f_face_surf = mq->i_f_face_surf;

    for (cs_lnum_t f_id = 0; f_id < n_i_faces; f_id++) {
      cs_real_t f_porosity = i_face_porosity[f_id];
      for (cs_lnum_t i = 0; i < 3; i++)
        i_f_face_normal[f_id][i] = f_porosity * i_face_normal[f_id][i];
      i_f_face_surf[f_id] = cs_math_3_norm(i_f_face_normal[f_id]);
    }

  }

  if (mq->b_f_face_surf != mq->b_face_surf) {

    const cs_real_3_t *restrict b_face_normal
      = (const cs_real_3_t *restrict)mq->b_face_normal;
    cs_real_3_t *restrict b_f_face_normal
      = (cs_real_3_t *restrict)mq->b_f_face_normal;
    cs_real_t *b_f_face_surf = mq->b_f_face_surf;

    for (cs_lnum_t f_id = 0; f_id < n_b_faces; f_id++) {
      cs_real_t f_porosity = b_face_porosity[f_id];
      for (cs_lnum_t i = 0; i < 3; i++)
        b_f_face_normal[f_id][i] = f_porosity * b_face_normal[f_id][i];
      b_f_face_surf[f_id] = cs_math_3_norm(b_f_face_normal[f_id]);
    }

  }

  BFT_FREE(i_face_porosity);
  BFT_FREE(b_face_porosity);

  /* Four set face factor */

  if (mq->i_f_face_factor != NULL) {
    const cs_lnum_2_t *i_face_cells
      = (const cs_lnum_2_t *)m->i_face_cells;
    const cs_real_t *i_face_surf = mq->i_face_surf;
    const cs_real_t *i_f_face_surf = mq->i_f_face_surf;

    for (cs_lnum_t f_id = 0; f_id < m->n_i_faces; f_id++) {
      cs_lnum_t ii = i_face_cells[f_id][0];
      cs_lnum_t jj = i_face_cells[f_id][1];

      cs_real_t face_porosity
        = CS_MAX(i_f_face_surf[f_id] / i_face_surf[f_id],
                 cs_math_epzero);

      mq->i_f_face_factor[f_id][0] = cell_porosity[ii] / face_porosity;
      mq->i_f_face_factor[f_id][1] = cell_porosity[jj] / face_porosity;
    }
  }