cs_cdo_quantities.h

Go to the documentation of this file.

#ifndef CS_CDO_QUANTITIES_H
#define CS_CDO_QUANTITIES_H
 
/*============================================================================
 * Manage geometrical quantities needed in CDO/MAC schemes
 *============================================================================*/
 
/*
  This file is part of code_saturne, a general-purpose CFD tool.
 
  Copyright (C) 1998-2026 EDF S.A.
 
  This program is free software; you can redistribute it and/or modify it under
  the terms of the GNU General Public License as published by the Free Software
  Foundation; either version 2 of the License, or (at your option) any later
  version.
 
  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  details.
 
  You should have received a copy of the GNU General Public License along with
  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
  Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
 
/*----------------------------------------------------------------------------*/
 
/*----------------------------------------------------------------------------
 *  Local headers
 *----------------------------------------------------------------------------*/
 
#include "base/cs_base.h"
#include "base/cs_math.h"
#include "cdo/cs_cdo_connect.h"
#include "cdo/cs_flag.h"
#include "mesh/cs_mesh.h"
#include "mesh/cs_mesh_quantities.h"
 
/*============================================================================
 * Macro definitions
 *============================================================================*/
 
/*============================================================================
 * Type definitions
 *============================================================================*/
 
typedef enum {
 
  CS_CDO_QUANTITIES_MEANV_CENTER,
  CS_CDO_QUANTITIES_BARYC_CENTER,
  CS_CDO_QUANTITIES_SATURNE_CENTER,
  CS_CDO_QUANTITIES_SUBDIV_CENTER,
 
} cs_cdo_quantities_cell_center_algo_t;
 
typedef enum {
 
  /* Set of geometrical quantities related to CDO schemes */
 
  CS_CDO_QUANTITIES_EB_SCHEME  = 1 << 0, /* =   1 */
  CS_CDO_QUANTITIES_FB_SCHEME  = 1 << 1, /* =   2 */
  CS_CDO_QUANTITIES_HHO_SCHEME = 1 << 2, /* =   4 */
  CS_CDO_QUANTITIES_VB_SCHEME  = 1 << 3, /* =   8 */
  CS_CDO_QUANTITIES_VCB_SCHEME = 1 << 4, /* =  16 */
  CS_CDO_QUANTITIES_CB_SCHEME  = 1 << 5, /* =  32 */
  CS_CDO_QUANTITIES_MAC_SCHEME = 1 << 6, /* =  64 */
 
} cs_cdo_quantities_bit_t;
 
/* Structure storing information about variation of entities across the
   mesh for a given type of entity (cell, face and edge) */
 
typedef struct {
 
  /* Measure is either a volume for cells, a surface for faces or a length
     for edges */
 
  double   meas_min;  /* Min. value of the entity measure */
  double   meas_max;  /* Max. value of the entity measure */
  double   h_min;     /* Estimation of the min. value of the diameter */
  double   h_max;     /* Estimation of the max. value of the diameter */
 
} cs_quant_info_t;
 
/* For primal vector quantities (edge or face) */
 
typedef struct {
 
  double  meas;       /* length or area */
  double  unitv[3];   /* unitary vector: tangent or normal to the element */
  double  center[3];
 
} cs_quant_t;
 
typedef struct cs_cdo_quantities_t { /* Specific mesh quantities */
 
  /* Keep the information about the removal of boundary faces in case of 2D
     computations */
 
  bool             remove_boundary_faces;
 
  /* Global mesh quantities */
 
  double           vol_tot;
 
  /* Cell-based quantities */
  /* ===================== */
 
  cs_lnum_t         n_cells;        // Local number of cells
  cs_gnum_t         n_g_cells;      // Global number of cells
  cs_real_3_t      *cell_centers;   // May be shared according to options
  cs_real_t        *cell_vol;       // Can be shared with cs_mesh_quantities_t
 
  cs_quant_info_t   cell_info;
 
  /* Face-based quantities */
  /* ===================== */
 
  cs_lnum_t         n_faces;        /* n_i_faces + n_b_faces */
  cs_lnum_t         n_i_faces;      /* Local number of interior faces */
  cs_lnum_t         n_b_faces;      /* Local number of border faces */
  cs_gnum_t         n_g_faces;      /* Global number of faces */
 
  /* Remark: cs_quant_t structure attached to a face (interior or border) can
     be built on-the-fly calling the function cdoq.get_face(f_id).
     See \ref get_face for more details.
 
     In order to reduce the memory consumption one shares face quantities with
     the ones defined in the legacy part and stored in the cs_mesh_quantities_t
     structure that's why a distinction is made between interior and border
     faces.
 
     cs_nvec3_t structure associated to a face can also be built on-the-fly
     using cdoq.get_face_nvec(f_id).
     See \ref get_face_nvec for more details.
  */
 
  cs_nreal_3_t     *i_face_u_normal;  // Can be shared with cs_mesh_quantities_t
  cs_real_t        *i_face_normal;    // Can be shared with cs_mesh_quantities_t
  cs_real_t        *i_face_surf;      // Can be shared with cs_mesh_quantities_t
  cs_real_3_t      *i_face_center;    // Can be shared with cs_mesh_quantities_t
  const cs_real_t  *i_dist;           // Always shared with cs_mesh_quantities_t
 
  cs_nreal_3_t     *b_face_u_normal;  // Can be shared with cs_mesh_quantities_t
  cs_real_t        *b_face_normal;    // Can be shared with cs_mesh_quantities_t
  cs_real_t        *b_face_surf;      // Can be shared with cs_mesh_quantities_t
  cs_real_3_t      *b_face_center;    // Can be shared with cs_mesh_quantities_t
  const cs_real_t  *b_dist;           // Always shared with cs_mesh_quantities_t
 
  cs_flag_cartesian_axis_t *face_axis; /* Enum for normal direction of faces
                                        * Not always allocated
                                        */
 
  /* Remark: cs_nvec3_t structure attached to a dual edge can be built
     on-the-fly to access to its length and its unit tangential vector using
     the function cdoq.get_dedge_nvec(shift)
 
     One recalls that a dual edge is associated to a primal face and is shared
     with two cells for an interior face and shared with one cell for a
     boundary face.  Scan this quantity with the c2f connectivity.
  */
 
  cs_real_t        *dedge_vector;   /* Allocation to 3*c2f->idx[n_faces] */
 
  cs_real_t        *pvol_fc;        /* Portion of volume surrounding a face
                                     * in each cell. This is a pyramid of
                                     * base the face and apex the cell center
                                     * Scanned with the c2f adjacency.
                                     * Not always allocated.
                                     */
  cs_quant_info_t   face_info;
 
  /* Edge-based quantities */
  /* ===================== */
 
  cs_lnum_t         n_edges;        /* Local number of edges */
  cs_gnum_t         n_g_edges;      /* Global number of edges */
 
  cs_real_t        *edge_vector;    /* Allocation to 3*n_edges
                                       Norm of the vector is equal to the
                                       distance between two vertices.
                                       Unit vector is the tangential direction
                                       attached to the edge */
 
  /* For each edge e belonging to a cell c, the dual face is built from the
     contributions of two triangles s(x_c, x_f, x_e) and s(x_c, x_f', x_e) with
     the faces f and f' belonging to F_e \cap F_c
     Scan this quantity with the c2e connectivity */
 
  cs_real_t        *dface_normal;   /* Vector-valued normal for each dual face
                                     * inside a cell associated to an edge */
  cs_real_t        *pvol_ec;        /* Portion of volume surrounding an edge
                                     * in each cell. Scanned with the c2e
                                     * adjacency.
                                     * Not always allocated. */
 
  cs_quant_info_t   edge_info;
 
  /* Vertex-based quantities */
  /* ======================= */
 
  cs_lnum_t         n_vertices;      /* Local number of vertices */
  cs_gnum_t         n_g_vertices;    /* Global number of vertices */
 
  cs_real_t        *pvol_vc;         /* Part of the dual cell associated to a
                                      * vertex in each cell. These quantities
                                      * are scanned thanks to the c2v
                                      * adjacency structure */
 
  const cs_real_t  *vtx_coord;       /* Coordinates of the mesh vertices.
                                      * Always shared with the global cs_mesh_t
                                      * structure */
 
  /* Dual volume related to the dual cell associated in a one-to-one pairing to
   * each vertex. This quantity has been synchronized in case of parallel
   * computing. Size of the array = n_vertices. Not always allocated */
 
  cs_real_t        *dual_vol;
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  inline cs_real_t
  get_face_surf(const cs_lnum_t f_id) const
  {
    return (f_id < this->n_i_faces) ? this->i_face_surf[f_id]
                                    : this->b_face_surf[f_id - this->n_i_faces];
  };
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  inline cs_real_t
  get_bface_surf(const cs_lnum_t bf_id) const
  {
    assert(bf_id >= 0 && bf_id < this->n_b_faces);
    return this->b_face_surf[bf_id];
  };
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  inline const cs_nreal_t *
  get_bface_normal(cs_lnum_t bf_id) const
  {
    assert(bf_id >= 0 && bf_id < this->n_b_faces);
    return this->b_face_u_normal[bf_id];
  };
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  inline const cs_real_t *
  get_face_vector_area(cs_lnum_t f_id) const
  {
    if (f_id < this->n_i_faces) /* Interior face */
      return this->i_face_normal + 3 * f_id;
    else /* Border face */
      return this->b_face_normal + 3 * (f_id - this->n_i_faces);
  };
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  inline const cs_real_t *
  get_face_center(cs_lnum_t f_id) const
  {
    if (f_id < this->n_i_faces) /* Interior face */
      return this->i_face_center[f_id];
    else /* Border face */
      return this->b_face_center[f_id - this->n_i_faces];
  }
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  cs_quant_t
  get_face(cs_lnum_t f_id) const;
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  cs_nvec3_t
  get_face_nvec(cs_lnum_t f_id) const;
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  cs_quant_t
  get_edge_center(cs_lnum_t e_id, const cs_cdo_connect_t *topo) const;
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  cs_nvec3_t
  get_edge_nvec(cs_lnum_t e_id) const;
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  cs_nvec3_t
  get_dedge_nvec(cs_lnum_t shift) const;
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  void
  log_summary() const;
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  void
  dump() const;
 
  /*----------------------------------------------------------------------------*/
  /*----------------------------------------------------------------------------*/
 
  void
  check(int                         verb,
        const cs_mesh_quantities_t *mq,
        const cs_cdo_connect_t     *topo) const;
 
} cs_cdo_quantities_t;
 
/*============================================================================
 * Global variables
 *============================================================================*/
 
/*============================================================================
 * Public function prototypes
 *============================================================================*/
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
static inline double
cs_compute_area_from_quant(const cs_quant_t    qa,
                           const cs_real_t    *xb)
{
  const double  xab[3] = {xb[0] - qa.center[0],
                          xb[1] - qa.center[1],
                          xb[2] - qa.center[2]};
  const double  cp[3] = {qa.unitv[1]*xab[2] - qa.unitv[2]*xab[1],
                         qa.unitv[2]*xab[0] - qa.unitv[0]*xab[2],
                         qa.unitv[0]*xab[1] - qa.unitv[1]*xab[0]};
 
  return 0.5 * qa.meas * cs_math_3_norm(cp);
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
long long
cs_cdo_quantities_get_time_perfo(void);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_cdo_quantities_set(cs_flag_t   option_flag);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_cdo_quantities_set_algo_ccenter(cs_cdo_quantities_cell_center_algo_t   algo);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
cs_cdo_quantities_t *
cs_cdo_quantities_build(const cs_mesh_t             *m,
                        const cs_mesh_quantities_t  *mq,
                        const cs_cdo_connect_t      *topo);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
cs_cdo_quantities_t *
cs_cdo_quantities_free(cs_cdo_quantities_t   *cdoq);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_cdo_quantities_free_cell_vol(cs_cdo_quantities_t *cdoq);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
const cs_real_t *
cs_cdo_quantities_get_pvol_fc(const cs_cdo_quantities_t *cdoq,
                              const cs_adjacency_t      *c2f);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_cdo_quantities_compute_pvol_fc(const cs_cdo_quantities_t *cdoq,
                                  const cs_adjacency_t      *c2f,
                                  cs_real_t                **p_pvol_fc);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
const cs_real_t *
cs_cdo_quantities_get_pvol_ec(const cs_cdo_quantities_t *cdoq,
                              const cs_adjacency_t      *c2e);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_cdo_quantities_compute_pvol_ec(const cs_cdo_quantities_t   *cdoq,
                                  const cs_adjacency_t        *c2e,
                                  cs_real_t                  **p_pvol_ec);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
const cs_real_t *
cs_cdo_quantities_get_dual_volumes(cs_cdo_quantities_t      *cdoq,
                                   const cs_cdo_connect_t   *connect);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_cdo_quantities_compute_dual_volumes(const cs_cdo_quantities_t   *cdoq,
                                       const cs_cdo_connect_t      *connect,
                                       cs_real_t                  **p_dual_vol);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_cdo_quantities_compute_i_tef(const cs_cdo_connect_t       *connect,
                                const cs_cdo_quantities_t    *cdoq,
                                cs_lnum_t                     f_id,
                                cs_real_t                     tef[]);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_cdo_quantities_compute_b_tef(const cs_cdo_connect_t       *connect,
                                const cs_cdo_quantities_t    *cdoq,
                                cs_lnum_t                     bf_id,
                                cs_real_t                     tef[]);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_cdo_quantities_compute_i_wvf(const cs_cdo_connect_t       *connect,
                                const cs_cdo_quantities_t    *cdoq,
                                cs_lnum_t                     f_id,
                                cs_real_t                     wvf[]);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_cdo_quantities_compute_b_wvf(const cs_cdo_connect_t       *connect,
                                const cs_cdo_quantities_t    *cdoq,
                                cs_lnum_t                     bf_id,
                                cs_real_t                     wvf[]);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_quant_dump(FILE *f, cs_lnum_t num, const cs_quant_t q);
 
#endif /* CS_CDO_QUANTITIES_H */